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I.  INTRODUCTION 


The  need  for  better  and  more  efficient  gas  turbines  requires  the  availability  of  cheap 
and  reliable  design  tools  for  blades  used  in  compressors  and  turbines.  Computational 
methods  are  the  preferred  choice  for  such  a  design  tool,  considering  the  cost  and  com¬ 
plexity  of  wind  tunnel  experiments. 

Among  the  computational  methods  available  today,  the  logical  choice  seems  to  be 
a  computer  code  that  can  solve  directly  the  full  Xavier  Stokes  equations.  However,  given 
the  state  of  the  art  in  both  algorithms  and  computer  hardware,  such  Xavier  Stokes 
solvers  are  restricted  only  to  supercomputers,  and  even  then  the  computation  time  is 
quite  long. 

In  order  to  enable  fast  and  efficient  computations,  the  viscous  in  viscid  interaction 
code  was  developed  by  Cebeci  [Ref.  1].  The  approach  used  in  this  code  is  to  solve  the 
outer  flow  field  using  potential  methods,  and  solving  the  boundary  layer  flow  using  a 
boundary  layer  method  subject  to  an  interaction  law,  that  couples  the  inner  and  the 
outer  flows.  This  interaction  law  is  needed  because  classical  boundary  layer  methods  fail 
in  areas  of  flow  reversal  and  separation,  which  are  very  common  in  real  life  flows. 

The  viscous  inviscid  interaction  code  was  originally  developed,  and  successfully  used, 
for  (lows  about  single  airfoils.  It  was  later  adapted  to  cascade  flows. 

In  this  thesis  the  applicability  of  the  code  to  cascades  was  investigated  by  comparing 
its  results  to  experimental  data.  It  was  found  that  although  the  code  can  reasonably 
predict  experimental  results  in  some  cases,  it  still  needs  improvements  before  it  can  be 
applied  generally  as  a  reliable  design  tool. 

A  major  restriction  in  improving  the  code  is  the  lack  of  a  wide  data  base  of  appro¬ 
priate  experimental  results.  Some  of  the  key  elements  in  the  code,  like  transition  and 
turbulence  modelling,  are  based  on  empirical  correlations,  and  more  detailed  and  accu¬ 
rate  experiments  should  be  performed,  before  a  better  understanding  of  these  phenom¬ 
ena  can  be  achieved. 

In  the  following,  the  theoretical  background  of  the  code  is  presented  in  Chapter  II, 
a  description  of  the  code  in  Chapter  III,  the  results  and  discussions  are  presented  in 
Chapter  IV  and  the  conclusions  and  recommendations  in  Chapter  V.  A  listing  of  the 
computer  program  is  given  in  Appendix  A. 
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II.  CASCADE  FLOW  PROBLEM  FORMULATION 


This  chapter  outlines  the  theoretical  background  of  the  viscous  in  viscid  interactive 
method  used  in  the  computer  code  to  investigate  cascade  flows.  Only  the  major  steps  in 
the  mathematical  developments  will  be  described  here.  A  detailed  description  of  the 
theory  and  the  numerical  methods  is  given  by  Cebeci  and  Bradshaw  [Ref.  2]  and  by 
Krainer  [Ref.  3[  on  which  this  chapter  is  based. 

A.  INYISCID  FLOW  METHOD 

Inviscid  flow  is  the  first  building  block  of  the  flow  and  is  solved  using  the  panel 
method.  The  incompressible  two  dimensional  outer  flow  must  satisfy  the  Laplace 
equation: 

Y24>  =  0 , 

subject  to  the  boundary  conditions  on  the  surface  of  the  blade: 


where  the  commonly  used  boundary  condition  of  zero  normal  velocity  on  the  surface  is 
replaced  by  a  specified  blowing  velocity  v„  to  allow  for  the  effect  of  the  boundary  layer 
on  the  outer  flow. 

In  addition,  the  Kutta  condition  must  be  satisfied,  in  order  to  prevent  the  existence 
of  discontinuous  pressure  distribution  near  the  trailing  edge  of  the  blade. 

Since  the  Laplace  equation  is  linear,  a  solution  to  a  complex  flow  field  can  be  built 
by  superposition  of  solutions  of  elementary  flows.  The  elementary  flows  used  in  the 
panel  method  are  the  uniform  parallel  flow  and  flows  about  singularities  (sources  and 
vortices). 

The  panel  method  is  based  on  replacing  each  blade  by  a  distribution  of  sources  and 
vortices  on  its  surface.  The  surface  is  divided  into  a  finite  number  of  straight  segments, 
called  panels. 

On  each  panel,  a  uniform  source  distribution  and  a  uniform  vorticity  distribution  is 
assumed.  The  source  strength  at  each  panel  is  set  to  satisfy  the  boundary  condition  at 
the  midpoint  of  the  panel  (called  the  control  point),  and  so,  in  general  the  source 
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strength  will  \ary  from  panel  to  panel.  The  vorticity  strength  is  assumed  to  he  the  same 
for  all  the  panels  and  is  set  to  satisfy  the  Kutta  condition. 

The  cascade  is  defined  as  an  infinite  row  of  similar  blades,  each  one  modelled  bv 
panels  of  source  and  vortex  distributions.  The  flow  at  each  point  is  found  by  summing 
the  contributions  of  all  the  singularities  on  all  the  blades  .  and  the  uniform  parallel  flow. 

A  useful  concept  in  dealing  with  such  flows  is  the  concept  of  influence  coefficients. 
\n  influence  coefficient  is  defined  as  the  velocity  at  a  point  induced  by  a  unit  strength 
singularity  placed  at  some  other  field  point.  The  influence  coefficients  are  a  function  of 
geometry  and  so  can  be  computed  for  a  given  cascade  and  a  given  choice  of  panel  ge¬ 
ometry. 

Using  the  influence  coefficients,  the  normal  and  tangential  velocities  at  each  control 
point  can  be  written  as  a  function  of  the  unknown  source  strength  of  each  panel  and  the 
unknown  vortex  strength.  Using  the  boundary  conditions,  by  equating  the  normal  ve¬ 
locity  at  each  control  point  to  the  prescribed  blowing  velocity  v„,  and  using  the  Kutta 
condition  (which  requires  equal  velocities  on  the  upper  and  on  the  lower  panels  at  the 
trailing  edge),  a  system  of  linear  equations  is  constructed. 

By  solving  the  system  of  equations,  the  strength  of  the  sources  and  vortices  is  found, 
and  the  velocities  (and  the  pressures)  can  be  computed  everywhere  in  the  flow  field. 

The  velocity  distribution  on  the  surface  of  the  blade,  computed  by  the  panel  method, 
is  used  as  the  boundary  condition  for  the  boundary  layer  flow  calculations. 

It  should  be  noted  that  in  the  panel  method  as  used  in  the  present  computer  code, 
there  is  no  modelling  of  the  wake,  and  its  effect  on  the  flow  field  is  ignored. 

B.  VISCOUS  FLOW  METHOD 

Viscous  flow  is  the  second  building  block  of  the  cascade  flow  and  it  is  applied  to  the 
thin  boundary  layer  near  the  blade  surface. 

1.  Boundary  Layer  Theory 

The  boundary  layer  concept,  first  suggested  by  L.  Prandtl,  assumes  that  the  flow 
field  can  be  divided  into  an  outer  flow  where  the  viscous  effects  are  negligible  compared 
to  inertia  effects,  and  a  thin  layer  close  to  the  surface  where  the  viscous  effects  cannot 
be  neglected.  The  complete  presentation  of  the  boundary  layer  theory,  and  the  devel¬ 
opment  of  the  boundary  layer  equations,  is  given  by  Schlichting  [Ref.  4). 

Under  the  assumptions  of  two  dimensional  incompressible  thin  boundary  layer 
flow,  the  Navicr  Stokes  equations  and  the  continuity  equation  reduce  to: 
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V 


cu 

ex 


+ 


cv 

Cy 

due 

dx 


=  0, 


with  the  boundary  conditions: 

.v  =  0  u(a,0)  =  0  ,  v(.t,0)  =  0  , 

y=ye  u(xye)  =  ue(x), 

V. 

where  b  denotes  1  + 


Writing  the  velocity  components  in  terms  of  a  stream  function  Y  : 


6\p 


_ 8x l>_ 

cx 

This  eliminates  the  continuity  equation  {which  the  stream  function  satisfies  by  defi¬ 
nition). 


Introducing  the  Falkner  Skan  transformation: 


V  =  V 


/  t*e 
VX 


A*y)  = 


l 


Juevx 


'Pixy) , 


the  momentum  equation  and  the  boundary'  conditions  transform  to: 


(bf"Y  +  +  m[l  -  CT')2]  =  x(f  -f^- -r  ~  , , 


>7  =  0  /'(jf,0)  =  0,J[xy)  =  0, 
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V  =  fie 


f'{x,  rje)  =  1  , 


where  m  is  defined  by: 

x  due 

m  =  ~e  —  • 

The  third  order  differential  equation  can  be  reduced  to  a  system  of  first  order 
differential  equations  by  introduction  of  two  new  variables  U  and  V  : 

L'  =  /' , 

V  =  L"  , 

(b\J  +  +  m(  1  -  L'2)  = 

2  \  cx  cx  ) 

with  the  boundary  conditions: 

V  =  0  L'(x,0)  =  0,J{x,0)  =  0, 

V  =  Ve  L’(^.  Ve)  -  1  • 

The  next  step  is  to  use  a  finite  difference  approach  to  solve  the  equations.  The 
box  method  is  applied  using  central  differencing  in  both  the  x  and  rj  directions,  and 
satisfying  the  equations  midway  between  nodes. 

Applying  the  box  method  results  in  a  system  of  nonlinear  equations  in  the  un¬ 
known  variables  (which  are  f,  U  and  V  in  each  node  along  the  rj  direction  at  the  current 
x  station). 

In  order  to  solve  the  nonlinear  system  the  Newton  iterative  procedure  is  used, 
linearizing  the  equations  first  about  the  solution  at  the  adjacent  upstream  station,  and 
then  about  the  preceding  iteration.  The  linearization  is  performed  by  letting: 

fr =fr'+sf;K, 

uj ,K  =  Uj’K_1  +  d\jjK  , 

VljK  =  WljK-]  +SV1JK  , 
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where: 


•  i  denotes  location  in  the  x  direction 

•  j  denotes  location  in  the  y  (  >/  )  direction 

•  k  indicates  the  iteration  counter 

This  linearization  results  in  a  system  of  linear  equations  for  the  unknown  in¬ 
crements:  3/}'w ,  and  3V‘-’ . 

This  system  of  equations  is  solved  repeatedly  until  the  changes  in  the  unknowns 
are  small  enough.  Since  the  system  is  block  tridiagonal,  Keller's  block  elimination 
method  is  used. 

The  method  described  so  far.  is  a  direct  boundary  layer  method.  It  can  be  used 
as  long  as  the  flow  does  not  separate.  Whenever  separation  or  flow  reversal  occurs,  and 
a  zero  skin  friction  coefficient  is  encountered,  the  equations  become  singular  and  the 
calculations  will  break  down. 

2.  Interactive  Boundary  Layer  Method 

The  interactive  boundary  layer  method  is  designed  to  overcome  the  difficulties 
encountered  at  regions  of  flow  reversal  and  separations.  In  such  areas  the  external  ve¬ 
locity  is  substantially  changed  by  the  viscous  effects  and  can  no  longer  be  considered  as 
a  known  boundary  condition  for  the  boundary  layer  flow. 

The  general  approach  to  the  solution  is  the  same  as  for  the  direct  method  but, 
since  the  outer  flow  is  unknown,  the  velocity  at  the  edge  of  the  boundary  layer  is  written 
as: 


ve)  = 


4e\ 


+  ir 


d  ,  ,V  dl 


where: 

1.  u,(x,y,)  is  the  total  velocity  at  the  edge  of  the  boundary  layer. 

2.  u,,(jc)  is  the  velocity  as  computed  by  the  inviscid  method. 

r 


3. 


1 

7T 


d\ 

x-l 


is  the  Hilbert  integral. 
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The  numerical  solution  of  the  boundary  layer  equations  follows  the  same  steps 
as  for  the  direct  method,  but  with  some  changes. 

The  transformation  of  the  stream  function  and  they?  coordinate  uses  a  constant 
velocity  ^  as  a  scaling  factor,  and  a  scaled  velocity  w  is  introduced: 


/  “o 

*  -  v  "v3r  y  ’ 


Ax,  v)  = 


1 


s/u0vx 


'J'(xj) , 


w  = 


ue(x  *>') 


Using  this  transformation,  the  boundary  layer  equations  become  a  system  of 
first  order  differential  equations: 

-  f  =  U, 

U'  =  V, 

W'  =  0, 


(bvy  +  -}rjv  +  x w4~  =  x(u4r--V-¥-  }, 

l  cx  \  cx  cx 


with  the  boundarv  conditions: 


V  =  0  U(*,0)  =  0  ,./(*, 0)  =  0  , 


V  =  Ve 


U(x,  Ve)  =  W(*>  Ve)  • 


W(x,  t]e)  = 


“o 


,  JL  [  J_(  [Kr^tr  s  dl 

n  }  (  V  «0  ^ ’  V e)V e  A£<Ve)^J  x—£  ’ 


The  finite  difference  box  method  is  used  to  solve  the  equations,  in  the  same  way 
as  it  was  used  for  the  direct  case,  but  with  two  additions: 
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1.  In  areas  of  flow  reversal  the  term  ucu'cx  is  omitted  to  assure  stable  integration 
(the  FLARE  approximation). 

2.  The  edge  velocity,  Wj  (where  J  denotes  the  edge  station)  which  involves  integration, 
is  approximated  by  : 

Wi  =  gi  +  c,7(VVj  >i}  —f5)  , 

where  g ...  and  c„  are  obtained  from  the  numerical  approximation  to  the  Hilbert  in¬ 
tegral  (which  will  be  presented  in  the  next  section  ). 

By  using  central  diffrencing  to  approximate  the  differential  equations,  a  svstem 
of  nonlinear  algebraic  equations  is  obtained  for  the  unknown  variables  (which  are 
f, ,  Uj,Vj  and  Wj).  To  solve  the  system  of  equations,  the  system  is  linearized  by  the 
Newton  iterative  procedure,  and  the  resulting  linear  system  is  solved  (for  the  new  un¬ 
known  variables  which  are  the  increments  bj ;-r ,  5 L')ir ,  8V';"  and  S W'-r ). 

The  solution  of  the  system  is  repeated  until  the  change  in  the  increments  is 
negligible  compared  to  the  preceding  iteration,  and  the  whole  process  is  performed  again 
at  the  next  downstream  station. 

3.  Interactive  Model 

The  interactive  model  is  used  to  couple  the  boundary  layer  to  the  external  flow. 
It  is  needed  in  areas  where  strong  interaction  occurs,  and  both  the  boundary  layer  and 
the  outer  flow  must  be  solved  simultaneously.  The  interaction  model  provides  the  outer 
boundary'  condition  to  the  boundary  layer  calculations  by  adding  a  correction  term  to 
the  external  velocity  computed  by  the  inviscid  flow'  method. 

The  external  velocity  is  assumed  to  consist  of  a  potential  flow’  term  (  u,,(x) )  and 
a  correction  term  due  to  viscous  effects  (  uti  (x) ): 


=  W<rlW  +  uei(x)  • 


The  viscous  effect  is  obtained  by  a  surface  distribution  of  sources  on  the  blade  (a  con¬ 
cept  first  suggested  by  Lighthill  [Ref.  5]).  The  normal  velocities  at  the  surface  of  the 
blade,  induced  by  these  sources,  displace  the  streamlines  from  the  surface  in  the  same 
way  that  the  actual  boundary  layer  displaces  them: 

db  (or)  v(x,  3  ) 
dx  ue{x)  ’ 

Where  v(x,  S’)  is  the  normal  velocity  at  the  displaced  surface. 
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Assuming  that  the  surface  can  be  approximated  to  be  a  Hat  plate,  the  normal 
velocity  will  be  half  the  local  source  strength  a  (x).  Assuming  also  that  the  inviscid  ve¬ 
locity  does  not  change  across  the  boundary  layer,  the  local  source  strength  will  be: 


2 


v(jc,0)  =  v(jc,  <5  )  — 


The  local  horizontal  velocity  induced  by  the  source  distribution,  is  the  cor¬ 
rection  term  to  the  inviscid  velocity,  and  can  be  represented  by  the  Hilbert  integral: 


x 

71 


rr(c) 


dc  =  ~~ 


rx. 


.X  —  C 


The  integration  is  carried  out  on  all  the  sources  on  the  surface,  since  the  horizontal 
velocity  is  influenced  by  all  the  sources. 

The  Hilbert  integral  is  then  approximated  by  a  finite  series: 


1 

7T 


Where  c,k  is  a  matrix  of  interaction  coefficients  which  arc  functions  of  the  geometry  only 
(  /  denotes  the  chordwise  position  where  uei  is  evaluated  and  k  is  the  location  of  the 
source  which  effects  ue6 ). 

Since  the  computation  of  uti  involves  values  of  <5*  downstream  of  the  current  x 
location,  which  are  not  known  yet,  these  terms  are  taken  from  the  previous  iteration 
using  a  relaxation  formula. 

4.  Turbulence  Model 

The  turbulence  model  used  here  is  the  algebraic  eddy  viscosity  formulation  of 
Cebeci  and  Smith  [Ref.  6].  According  to  the  model  used  in  the  present  computer  code, 
the  eddy  viscosity  v,  is  defined  by  two  different  expressions,  for  the  inner  region  and  for 
the  outer  region: 
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for  0  <>•  <  jc , 


v 


for  yc<y<5. 


Where: 


A  = 


26v 

v-^V/2 

cy  J 


l 


1  *f 

5.5(k/<5)6 

» 

0.0168 

1  -  J? 

du’dx  "1 

2.S 

CU  CV  ; 

a  = 


_ 6 _ 

1  +  2Rt(2  -  Rt) 


for  Rt  <  1 , 


1  4-  Rt 
Rt 


for  RT  >  1  , 


R-r 


(  -  u'V) 


max 


The  distance  from  the  wall  to  the  point  between  the  two  regions,  yc,  is  chosen  such  that 
the  viscosity  will  be  continuous. 
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The  intermittency  factor,  y„  is  defined  by: 


'hr  =  1  -  exp 


u. 


G/ 


~  x„) 


di_ 


Where: 

•  R,„  is  the  Reynolds  number  based  on  external  velocity  and  transition  location. 

•  Gr  is  an  empirical  constant,  originally  assigned  the  value  1200. 

Cebeci  and  Bradshaw  [Ref.  2,  p.246]  described  a  different  expression  for  the 
variable  A  in  the  inner  region  viscosity  formula: 


A  = 


26v 


(1  -  1I.8/T)1'2 )1/2 

°y  J  max 


Where: 

*  dl,e 

P  (v4Ljn  dx 

This  version  of  the  turbulence  model  was  not  implemented  in  the  original  computer 
code.  During  the  work  on  this  thesis,  the  effect  of  the  modified  turbulence  model  was 
investigated. 

A  different  intermittency  distribution  was  implemented  successfully  by  Rodi  and 
Schonung  [Ref.  7  [  for  transition  over  separation  bubbles.  They  used  for  Gr  the  ex¬ 
pression: 

c  =  100 

”  exp(0.997u)  ‘ 

Where  Tu  is  the  turbulence  level  in  the  free  flow.  This  intermittency  model  was  also  in¬ 
vestigated  during  the  work  on  this  thesis. 


5.  Transition 


The  prediction  of  transition  from  laminar  to  turbulent  flow  is  very  difficult  and 
has  to  rely  on  empirical  correlations.  The  relation  used  here  to  predict  the  onset  of 
transition  is  a  combination  of  Michel's  method  and  the  e9  method,  and  is  given  by 
Cebeci  and  Bradshaw  [Ref.  2  ,  p.  153]: 


v  = 


1.174  1  + 


22400  \  r»  0.46 


R* 


R. 


Where: 

1.  R„  is  the  Reynolds  number  based  on  the  momentum  thickness  at  the  onset  of 
transition. 

2.  R,  is  the  Revnolds  number  based  on  x  at  the  onset  of  transition. 

€Xtr  " 

In  the  computer  code,  if  a  laminar  separation  is  detected  before  transition  oc¬ 
curs,  the  onset  of  transition  is  assumed  at  the  point  of  laminar  separation. 
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III.  DESCRIPTION  OF  THE  COMPUTER  CODE 


The  computer  code  used  here  to  investigate  cascade  flows  was  written  by  Cebeci, 
and  is  based  on  the  numerical  formulation  that  was  outlined  in  the  previous  chapters. 
In  this  chapter  the  general  structure  and  the  major  subroutines  of  the  code  will  be  de¬ 
scribed. 

A.  GENERAL  STRUCTURE  OF  THE  MAIN  PROGRAM 

The  main  program  reads  in  the  cascade  data  (blade  coordinates,  spacing  and  stagger 
angle),  the  flow  data  (inlet  angle  and  Reynolds  number),  and  transition  parameters.  The 
transition  onset  on  each  surface  of  the  blade  can  be  computed  by  the  program,  or  can 
be  input  by  the  user.  The  intermittency  parameter  G  should  be  specified  by  the  user. 

The  program  then  calls  subroutine  POTNL  to  compute  the  outer  inviscid  flow  field 
for  the  first  cycle.  The  output  of  subroutine  POTXL  is  the  external  velocity  distribution 
on  the  surface  of  the  blades.  This  velocity  distribution  is  then  transferred  to  subroutine 
CASBLP,  which  calculate  the  boundary  layer  flow. 

Subroutine  CASBLP  returns  the  displacement  thickness  distribution  and  the  blow¬ 
ing  velocity  distribution  on  the  blades  to  the  main  program.  This  data  is  then  transferred 
back  to  subroutine  POTNL  to  the  next  cycle  of  calculations. 

The  program  repeats  the  cycles  of  calculations  by  calling  the  two  subroutines,  until 
the  specified  number  of  cycles  is  reached,  or  until  a  convergence  criterion  is  satisfied. 

B.  DESCRIPTION  OF  THE  SUBROUTINES 
1.  Subroutine  POTNL 

This  subroutine  solves  the  inviscid  outer  flow  by  using  the  panel  method.  The 
subroutine  calculates  the  influence  coefficients  and  calculates  the  velocities  subject  to  the 
boundary  conditions. 

The  velocities  are  evaluated  on  the  displaced  surface  (the  surface  created  by 
adding  the  displacement  thickness  to  the  original  surface  of  the  blade).  The  input  to  this 
subroutine  includes  the  cascade  geometry,  the  blowing  velocity  and  the  displacement 
thickness  (for  the  first  cycle  both  the  displacement  thickness  and  the  blowing  velocity 
are  taken  to  be  zero). 
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2.  Subroutine  CASBLP 

This  subroutine,  called  by  the  MAIN  program,  receives  the  blade  geometry  and 
the  velocity  distribution  as  input. 

It  transforms  the  x.y  blade  coordinates  to  the  chordwise  tangential  coordinates 
and  smooths  the  velocity  data  (during  the  work  on  this  thesis  it  was  found  that 
smoothing  the  velocity  data  prevents  the  detection  of  the  separation  bubble  near  the 
leading  edge,  and  therefore  it  was  eliminated).  The  subroutine  then  calls  subroutine 
COMPBL  for  further  calculations. 

3.  Subroutine  COMPBL 

This  subroutine  finds  the  stagnation  point  and  controls  the  generation  of  the 
boundary  layer  calculation  grid  for  each  surface  (the  grid  starts  at  the  stagnation  point 
and  includes  91  points  in  the  chordwise  direction  for  the  upper  surface  and  71  points  on 
the  lower  surface). 

The  subroutine  then  calls  subroutine  BL2D  which  calculates  the  boundary  layer 
parameters  for  each  surface  (BL2D  is  called  twice,  first  for  the  upper  surface  and  then 
for  the  lower  surface). 

4.  Subroutine  BL2D 

This  subroutine  computes  the  displacement  thickness  and  the  blowing  velocity 
and  returns  them  back  to  the  calling  subroutine  (COMBL)  in  arrays  compatible  with  the 
potential  flow  calculations  (one  array  that  contains  all  the  points  of  the  blade,  first  the 
lower  surface  starting  at  the  trailing  edge  and  proceeding  forward,  and  then  the  upper 
surface,  starting  at  the  leading  edge  and  proceeding  backwards). 

BL2D  calls  the  following  subroutines: 

1.  Subroutine  INPUT  which  calculates  the  following: 

a.  NS,  the  switching  point  between  direct  and  interactive  boundary  layer  calcu¬ 
lations  (this  point  is  set  at  the  first  pressure  peak  when  the  blade  is  scanned  from 
leading  edge  towards  the  trailing  edge) 

b.  NTR,  transition  location  (only  if  the  transition  location  is  an  input.  Otherwise 
it  is  calculated  by  subroutine  TRNS). 

c.  GMTR,  the  distribution  of  the  intermittency  factor  y„. 

In  addition  this  subroutine  generates  the  boundary  layer  grid  in  the  rj  direction  and 

the  initial  velocity  profile,  by  calling  subroutine  INTL. 

2.  CALCIJ,  calculates  the  c„  coefficients  used  in  the  Hilbert  integral  approximation. 

3.  EDDY,  calculates  the  eddy  viscosity  (called  only  after  transition  has  been  de¬ 
tected). 
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4.  COEFTR,  calculates  the  coefficients  of  the  boundary  layer  finite  difference 
equations  in  transformed  form  (  for  the  direct  method  calculations). 

5.  SOLVE3,  solves  the  linearized  boundary  layer  equations  for  the  F,U  and  V  vari¬ 
ables  by  computing  the  increments  <5F,  dU  and  (5V 

The  subroutine  then  checks  the  convergence  of  the  Newton  iterations  and  re¬ 
peats  the  calculations  if  needed.  If  the  subroutine  detects  flow  separation  or  if  it  reaches 
the  switching  point  NS,  subroutine  MAIN2  is  called  for  the  interactive  method  calcu¬ 
lations.  Otherwise,  the  subroutine  proceeds  to  the  next  chordwise  point  of  the  grid  (NX) 
and  repeats  the  calculations. 

5.  Subroutine  MAIN2 

This  subroutine  calculates  the  boundary'  layer  parameters  by  the  interactive 
method.  The  subroutine  performs  the  following  steps: 

1.  It  first  calls  subroutines  JOINT  and  CO.MGI  to  compute  the  interaction  coeffi¬ 
cients. 

2.  In  regions  of  laminar  flow  it  calls  the  following  subroutines: 

a.  COEF,  which  calculates  the  coefficients  of  the  boundary  layer  finite  differences 
equations. 

b.  SOLV4,  solves  for  the  variables  F,  U,  V  and  W  by  computing  the  increments 
<5F  ,  dU  ,  <5V  and  d\V  . 

c.  TRANS,  to  check  if  the  condition  for  transition  is  satisfied  (it  also  checks  for 
laminar  separation  and  initiates  transition  at  the  point  of  laminar  separation  if 
it  is  detected). 

The  subroutine  then  checks  for  convergence  of  Newton  iterations  and  repeats  the 
calculations  as  needed. 

3.  In  regions  of  turbulent  flows  the  subroutine  calls  the  following  subroutines: 

a.  EDDY,  to  compute  the  eddy  viscosity  parameter  B  (B  =  1  +  vyv  ). 

b.  COEF  and  SOLV4,  the  same  as  for  laminar  flow. 


6.  Subroutine  OUTPUT 

This  subroutine  computes  the  boundary'  layer  parameters.  It  is  called  with  a 

parameter  "INDEX"  which  determines  the  type  of  calculations: 

1.  For  INDEX  =1  the  computations  relates  to  transformed  coordinates  (direct 
boundary  layer  method)  using  the  relations: 


c,  = 


2  V(1.2)  B(  1 ,2) 
/R 

N'  ‘Vx 
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Vw  =  nvfvC-2)- 
D  =  ueS  v  Re  , 


d  = 


0  = 


.t 


V  R* 


(rj(\P)-F(NP)), 


sp 


\  R-f 


E^c-u,) 


7=2 


Where  V(l,2)  and  B(  1 ,2 )  are  the  velocity  gradient  and  the  viscosity  parameter  at  the 
surface,  respectively,  and  j/(NP)  and  F(NP)  are  rj  and  F  evaluated  at  the  edge  of  the 
boundary  layer. 

2.  For  INDEX  =  2  the  subroutine  calculates  the  boundary  layer  parameters  for  semi- 
-transformed  coordinates  (interactive  boundary  layer  method)  using  the  relations: 

2  V(1.2)  B(  1 .2) 

Lf  ~  —  — 


V*  V  K  [W (SF)Y 
V(l,2) 


V„  - 


ue  =  U(A T) , 

D  =  (L  >?  -  F)v/7  , 

**  =  ("--fV7- 


For  NX  >  NTR  (after  transition  has  been  detected),  subroutine  SMPSON  is 
called  (subroutine  SMPSON  calculates  the  coefficient  of  the  outer  region  eddy  viscosity). 
The  subroutine  then  prints  out  the  velocity  profiles  at  the  required  stations. 

7.  Subroutine  TRANS 

This  subroutine  calculates  the  transition  location  based  on  the  Michel  criterion 
or  based  on  laminar  separation  (whichever  occurs  first).  If  transition  has  been  detected 
the  intermittency  distribution  is  calculated  for  all  the  remaining  points  of  the  surface. 

8.  Subroutine  FILLUP 

This  subroutine  increases  the  number  of  points  in  the  boundary  layer  grid  (in 
the  direction)  as  needed.  It  also  fills  up  the  arrays  of  F,  U,  B,  W  and  V  between  the 
edge  of  the  boundary  layer  to  the  end  of  the  arrays  (with  V  =  0,  W,B  and  U,  with  the  last 
values  that  they  had  in  the  edge  of  the  boundary  layer  and  F  as  the  integral  of  U). 
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9.  Subroutine  EDDY 

This  subroutine  calculates  the  eddy  viscosity  using  the  Cebeci--Smith  two  layer 
eddy  viscosity  formula.  It  receives  the  vectors  L,V  and  rj  at  a  point  and  computes  the 
viscosity  vector  B. 

10.  Subroutine  INTL 

This  subroutine  generates  the  boundary  layer  grid  in  the  r\  direction.  It  sets  the 
number  of  grid  points  and  generates  the  initial  velocity  profile. 
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IV.  RESULTS  AND  DISCUSSION 

The  viscous  inviscid  interaction  code  was  run  with  several  cascades  on  which  ex¬ 
perimental  data  is  available.  In  order  to  enable  a  thorough  comparison  between  exper¬ 
imental  results  and  the  computed  results,  a  very  detailed  experimental  data  base  is 
needed.  The  data  should  include  measurements  of  the  boundary  layer  development  along 
the  blade,  velocity  profiles  along  the  boundary  layer,  transition  location  and  distribution, 
flow  separation,  and  external  velocity  distribution. 

Unfortunately,  very  few  cascade  experiments  has  been  performed,  which  obtained 
the  required  data  with  sufficient  accuracy,  due  mostly  to  the  lack  of  appropriate  meas¬ 
urement  equipment.  Only  recently,  with  the  introduction  of  non-interfering  methods 
like  the  Laser  Doppler  Velocimeter  (LDV),  the  required  data  can  be  measured  accu¬ 
rately. 

Recently  an  experiment  involving  the  investigation  of  a  linear  compressor  cascade 
of  Controlled  Diffusion  Blading  (which  will  be  referred  here  as  the  CD  cascade)  has  been 
carried  out  by  Elazar  [Ref.  8],  Most  of  the  work  in  the  present  thesis,  involves  compar¬ 
ison  of  the  computer  code  results  with  Elazar's  experimental  results. 

Other  cascades  that  were  investigated  are: 

1.  A  shockless,  supercritical  airfoil  cascade,  designed  in  1974  by  Korn  in  cooperation 
with  Pratt  &  Whitney  Aircraft  (referred  here  as  the  P  &  W  cascade).  The  exper¬ 
imental  results  of  the  cascade  were  obtained  from  a  report  by  Hobbs,  Wagner. 
Dannenhoffer  and  Dring  [Ref.  9]. 

2.  Stator  blade  of  a  single  stage  axial  compressor  (referred  here  as  the  C4  cascade). 
The  blade  profile  is  the  British  C4  section  (10%  thickness)  on  a  circular  arc  camber 
line.  The  experiment  has  been  performed  by  Walker  [Ref.  10].  The  detailed 
boundary  layer  measurements  are  not  presented  in  the  report  and  were  obtained 
directly  from  the  author. 

The  code  failed  to  run  with  two  other  cascades: 

1.  A  highly  loaded,  double  circular  arc  blade  with  a  sharp  leading  edge  and  a  sharp 
trailing  edge,  used  in  a  compressor  cascade  that  was  investigated  by  Deutsch  and 
Zierk  [Ref.  11]. 

2.  V2  double  circular  arc  blade,  highly  loaded  cascade.  This  cascade  was  investigated 
by  Hoheisel  and  Seyb  [Ref.  12]. 

In  both  cases  the  code  calculated  the  potential  flow  successfully  but  failed  in  trying 
to  compute  the  first  cycle  of  the  boundary  layer  calculations. 
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A.  CD  CASCADE 

The  experimental  data  for  the  CD  cascade  was  obtained  at  M  =  0.25.  R,  =700000 
and  at  three  inlet  angles:  40°  (the  design  condition).  43.4°  and  46°  .  The  spacing  was 
0.6  of  the  chord  and  the  stagger  angle  14.27°  .  A  general  layout  of  the  cascade  is  shown 
in  Figure  1  on  page  20. 

The  following  observations  were  concluded  from  the  experiment: 

1.  A  separation  bubble  exists  near  the  leading  edge  on  the  upper  surface  at  all  the  inlet 
angles.  The  bubble  became  larger  at  increased  inlet  angles. 

2.  Transition  from  laminar  to  turbulent  flow  occurred  above  the  separation  bubble 
(on  the  upper  surface). 

3.  Transition  on  the  lower  surface  occurred  at  midchord. 

4.  The  boundary  layer  thickness  on  the  upper  surface  increased  with  inlet  angle,  and 
reached  a  thickness  of  15%  chord  at  the  highest  inlet  angle.  The  boundary  layer 
thickness  on  the  lower  surface  did  not  change  significantly  with  inlet  angle. 

5.  The  turbulent  boundary  layer  on  both  surfaces  remained  fully  attached  at  all  the 

inlet  angles. 

1.  Transition  location  and  intermittency  distribution 

The  effects  of  the  transition  location  and  the  intermittency  factor  were  investi¬ 
gated.  The  code  was  first  run  with  the  transition  location  calculated  by  the  code,  and 
with  several  values  of  the  intermittency  factor  Gr.  It  was  found  that  the  code  did  not 
run  with  G,  =  1200  (which  is  the  value  used  usually  for  high  Reynolds  numbers).  The 
highest  value  of  Gv  with  which  the  code  run  successfully  was  900. 

The  code  failed  to  predict  the  separation  bubble  on  the  upper  surface,  and  pre¬ 
dicted  laminar  separation  at  7S%  chord  on  the  lower  surface  (which  did  not  occur  in  the 
experiment).  Transition  on  the  upper  surface  occurred  at  41%  chord  (detected  by 
Michel's  criterion)  and  at  78%  chord  on  the  lower  surface  (at  laminar  separation). 

The  shape  factor  computed  by  the  code  was  compared  to  the  experimental  re¬ 
sults.  As  can  be  seen  in  Figure  2  on  page  21  the  shape  factor  as  predicted  by  the  code 
deviates  substantially  from  the  actual  results,  due  mainly  to  the  different  transition  lo¬ 
cation. 

On  the  lower  surface,  as  can  be  seen  in  Figure  3  on  page  22  the  shape  factor 
deviates  even  more  from  the  experimental  results.  In  this  figure  the  effect  of  changing 
the  intermittency  factor  G,  can  be  seen.  For  both  the  extreme  values  of  G, .  10  and  900, 
the  computed  shape  factor  curve  is  far  from  agreement  with  the  actual  results. 
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Figure  3.  Shape  factor  compai 
tiie  code  (/?  =  40°). 
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on  the  lower  surface:  Transition  computed  by 


Since  the  code  did  not  predict  the  existence  of'  the  separation  bubble  on  tire 
upper  surface,  it  was  decided  to  try  to  eliminate  the  smoothing  of  the  external  velocity 
as  computed  by  the  potential  flow  subroutine  (originally,  the  velocities  were  smoothed 
in  subroutine  CASBLP  prior  to  boundary  layer  calculations).  It  was  found  that  without 
smoothing  the  velocitic-'  a  small  separation  bubble  is  predicted  by  the  code  at  4%  of 
chord.  The  onset  of  transition  is  set  by  the  code  at  the  beginning  of  the  separation 
bubble. 

The  code  was  run  with  unsmoothed  velocities  with  two  values  of  Gy.  10  and  900. 
The  shape  factor  behavior  can  be  seen  in  Figure  4  on  page  24.  Changing  the  value  of 
Gv  did  not  change  the  shape  of  the  curve  much,  and  generally  the  shapes  of  the  com¬ 
puted  and  the  experimental  curves  look  alike. 

The  elimination  of  the  velocity  smoothing  in  the  code,  also  affects  the  thickness 
of  the  boundary  layer.  In  Figure  5  on  page  25  the  displacement  thickness  is  plotted  for 
both  cases  (with  and  without  the  velocity  smoothed).  Without  smoothing,  the  displace¬ 
ment  thickness  is  much  thicker,  especially  on  the  rear  half  of  the  blade,  which  is  closer 
to  the  actual  results. 

The  effect  of  changing  the  intermittency  distribution  to  the  one  used  by  Rodi 
and  Schonung  [Ref.  7]  was  investigated.  It  was  found,  as  can  be  seen  in  Figure  6  on  page 
26  that  the  effect  of  the  new  model  is  equivalent  to  using  G,  in  the  present  model. 

On  the  lower  surface  it  was  necessary  to  run  the  code  with  transition  as  input, 
to  get  reasonable  results,  as  can  be  seen  in  Figure  7  on  page  27  for  transition  input  at 
21° o  of  chord. 

At  the  off  design  conditions  (inlet  angles  of  43.4°  and  46°)  a  similar  behavior  of 
the  transition  has  been  observed,  as  can  be  seen  for  example  in  Figure  8  on  page  28  for 
the  upper  surface  and  in  Figure  9  on  page  29  for  the  lower  surface,  both  at  /?  =  46°. 


Figure  4.  Shape  factor  on  the  upper  surface  without  velocity  smoothing. 
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Figure  8.  Shape  factor  at  /J  =  46°  on  the  upper  surface 
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Figure  9.  Shape  factor  at  (l  =  46°  on  the  loner  surface. 


29 


2.  External  Velocity  Distribution 

The  external  velocity  distribution,  computed  by  the  code  using  the  interaction 
law,  was  compared  to  the  experimentally  measured  velocities.  It  was  found  that  in  gen¬ 
eral.  the  velocities  measured  experimentally,  were  higher  than  those  computed  by  the 
code  for  all  the  inlet  angles. 

There  are  two  possible  sources  to  the  discrepancy  in  the  velocities: 

1.  The  computer  code  calculates  pure  2-D  flows.  In  the  experiment  the  flow  was 
observed  to  accelerate  due  to  the  effect  of  the  boundary  layer  on  the  side  walls  (a 
3-D  effect).  This  effect  was  calculated  in  the  experiment  and  is  referred  to  as  the 
AYDR  correction  [Ref.  8,  p.43], 

2.  The  flow  accelerates  due  to  the  thickening  of  the  boundary  layer.  Since  the 
boundary’  layer  as  computed  by  the  code  is  substantially  thinner  than  the  actual 
boundary  layer  (as  will  be  discussed  in  the  next  section)  the  external  velocities 
predicted  by  the  code  arc  smaller. 

To  compensate  for  the  first  error  source,  all  the  computed  velocities  were  com¬ 
pared  to  the  experimental  velocities  corrected  by  the  AVDR  correction.  The  comparison 
between  the  velocities  can  be  seen  in  Figure  10  on  page  31  for  p  =  40°,  in  Figure  11 
on  page  32  for  /?  =  43.4°  and  in  Figure  12  on  page  33  for  /?  =  46°. 

It  can  be  seen  from  the  figures  that  the  difference  between  the  computed  and  the 
experimental  velocities  is  larger  on  the  lower  surface.  The  reason  might  be  the  method 
with  which  the  correction  to  the  inviscid  velocity  is  computed.  The  assumption  on  which 
the  interaction  law  is  based,  is  that  only  sources  (representing  the  viscous  effects)  on  the 
surface  being  considered,  affect  the  local  velocity.  In  reality,  the  boundary  layer  on  both 
surfaces  affects  the  local  velocity  (because  the  boundary  layer  developed  on  the  upper 
surface  of  a  blade,  causes  a  velocity  disturbance  that  is  felt  on  the  lower  surface  of  the 
adjacent  blade). 

Since  the  boundary  layer  on  the  lower  surface  is  much  thinner,  its  effect  on  the 
velocity  on  the  upper  surface  is  much  smaller  than  the  effect  of  the  upper  surface 
boundary  layer  on  the  lower  surface  velocity. 
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Figure  11.  External  velocity  at  ft  =  43.4° 
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3.  Boundary  Layer  Thickness 

The  boundary  layer  thickness  as  computed  by  the  code  was  compared  with  thr 
experimental  results  by  comparing  the  displacement  thicknesses. 

It  was  found  that  on  the  lower  surface  the  computed  and  the  actual  displace¬ 
ment  thickness  agree  quite  well,  as  can  be  seen  in  Figure  13  on  page  35  for  /?  =  40°, 
Figure  14  on  page  36  for  /?  =  43.4°  and  in  Figure  15  on  page  37  for  P  =  46°  . 

On  the  upper  surface  the  displacement  thicknesses  computed  by  the  program 
are  significantly  thinner  than  those  measured  experimentally.  The  difference  between  the 
computed  and  the  actual  thickness  increases  along  the  blade  and  it  increases  with  in¬ 
creased  inlet  angle.  It  was  found  that  by  using  a  different  expression  for  the  inner  region 
eddy  viscosity  (as  mentioned  in  chapter  II),  the  displacement  thickness  can  be  increased, 
but  the  difference  between  the  actual  and  the  computed  thickness  is  still  substantial,  es¬ 
pecially  at  the  higher  inlet  angles.  Figure  16  on  page  38,  Figure  17  on  page  39  and 
Figure  18  on  page  40  shows  the  displacement  thickness  on  the  upper  surface  for  the 
three  inlet  angles,  with  the  original  and  the  modified  eddy  viscosity  models. 

The  large  error  in  the  prediction  of  the  boundary  layer  thickness,  can  be  the  re¬ 
sult  of  several  reasons: 

1.  The  transition  model  used  in  the  code,  sets  the  onset  of  transition  at  the  first  point 
of  laminar  separation.  It  causes  rapid  transition  to  turbulent  flow  which  reattaches 
immediately,  resulting  in  a  very  small  separation  bubble  compared  to  the  bubble 
observed  in  the  experiment. 

2.  The  turbulent  model  used  in  the  code  could  be  inaccurate.  It  was  derived  based  on 
empirical  data  obtained  in  single  airfoil  experiments  and  not  with  cascades.  In  ad¬ 
dition  the  present  model  does  not  include  the  effects  of  the  free  stream  turbulence 
(that  was  relatively  high  in  the  experiment). 

3.  The  boundary  layer  as  measured  in  the  experiment  was  quite  thick,  especially  at  the 
higher  inlet  angles  (it  reached  15%  of  the  chord  at  P  =  46°).  Such  a  thick  bound¬ 
ary  layer  may  violate  the  basic  assumptions  on  which  the  boundary  layer 
equations,  and  the  interaction  law,  were  based  (especially  when  the  spacing  be¬ 
tween  the  blades  is  small,  60%  chord  in  this  case). 

It  was  suggested  that  one  of  the  possible  reasons  to  the  inaccurate  prediction 
of  the  boundary  layer  is  the  blunt  trailing  edge  of  the  blade,  that  might  cause  difficulties 
in  the  computations.  A  modified  blade,  with  a  sharp  trailing  edge  has  been  run,  and  the 
displacement  thickness  distribution  can  be  seen  in  Figure  19  on  page  41.  As  can  be  seen 
in  the  figure  the  sharp  trailing  edge  affects  only  the  boundary  layer  adjacent  to  the 
trailing  edge,  and  therefore  cannot  provide  an  explanation  to  the  difference  between  the 
actual  and  the  computed  displacement  thickness. 
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Figure  14.  Displacement  thickness  on  the  loner  surface  (/?  =  43.4°  ) 
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Figure  16.  Displacement  thickness  on  the  upper  surface  (/J  =  40°  ) 
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Figure  19.  The  effect  of  sharp  trailing  edge. 
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4.  Comparison  to  a  Navier  Stokes  Code 

A  limited  comparison  of  the  experimental  and  the  computed  results  with  a 
Navier  Stokes  (N.S.)  code  '’■’lculations  has  been  performed.  The  N.S.  code  has  been  de¬ 
veloped  and  run  by  S.  J.  Shamroth  of  Scientific  Research  Associates  inc.  in  cooperation 
with  Pratt  and  Whitney  Aircraft. 

Since  the  N.S.  code  does  not  compute  the  displacement  thickness,  the  velocity 
profiles  near  the  surface  of  the  blade  were  compared.  The  comparisons  were  made  at 
90°  o  chord  on  the  suction  surface  for  all  three  inlet  angles. 

At  the  design  point.  /?  =  40°,  shown  in  Figure  20  on  page  43.  both  the  inter¬ 
active  code  and  the  N.S.  code  failed  to  predict  accurately  the  actual  velocity  profile.  In 
this  case  the  interactive  code  seems  to  yield  somewhat  better  results  than  the  N.S.  code. 

At  the  higher  inlet  angles.  /?  =  43.4°  and  /?  =  46°,  shown  in  Figure  21  on  page 
44  and  in  Figure  22  on  page  45  respectively,  the  N.S.  calculations  show  significantly 
bettei  agreement  with  the  experimental  results  than  the  interactive  code. 

From  these  comparisons,  it  can  be  seen  that  the  interactive  code  deviation  from 
the  actual  results  increases  with  increased  inlet  angle  (increased  loading  of  the  cascade), 
whereas  the  N.S.  code  deviation  seems  to  decrease  with  increased  inlet  angle. 
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Figure  21.  The  results  of  the  N.  S.  code  at  P  =  43.4° 
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B.  P  &  W  CASCADE 

The  experimental  data  for  the  P  &.  W  cascade  was  obtained  at  inlet  flow  angle  of 
52°.  at  M  =  0.11.  and  Reynolds  number  of  478000.  The  cascade  had  a  stagger  angle  of 
15.75°  and  0.7  spacing.  A  general  layout  of  the  cascade  is  shown  in  Figure  23  on  page 
47. 

A  comparison  of  the  computed  and  the  measured  pressure  coefficients  on  the  blade 
is  shown  in  Figure  24  on  page  48.  There  is  a  good  agreement  between  the  computed 
and  the  measured  CP. 

The  displacement  thickness  was  measured  in  the  experiment  only  at  96.8%  of  chord. 
This  measurement  is  compared  to  the  computed  results  in  Figure  25  on  page  49.  As  can 
be  seen,  the  computed  and  the  measured  data  agree  almost  perfectly  on  the  lower  sur¬ 
face,  and  quite  well  on  the  upper  surface.  The  difference  observed  on  the  upper  surface 
is  caused  by  the  early  prediction,  by  the  code,  of  trailing  edge  separation,  a  short  dis¬ 
tance  upstream  of  the  actual  location.  This  can  also  be  observed  when  comparing  the 
velocity  profiles  at  that  point,  in  Figure  26  on  page  50.  The  computed  velocity  curve 
shows  a  small  zone  of  reversed  flow  near  the  surface  of  the  blade.  This  reversed  flow 
could  be  the  result  of  a  too  early  prediction  of  trailing  edge  separation  by  the  code,  or 
it  could  have  existed  in  the  actual  (low  but  not  detected  because  of  its  size. 
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Figure  23.  Pratt  &  Whiteny  cascade 
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Figure  26.  Velocity  profile  at  96.8%  chord  on  the  upper  surface. 
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C.  C4  CASCADE 


The  C4  cascade  has  a  stagger  angle  of  29.5°,  a  camber  angle  of  31.1°  and  spacing 
of  0. 992.  It  has  been  tested  at  Reynolds  numbers  of  about  200000,  and  inlet  angles  of 
34.  r  to  47.7°  (which  corresponds  to  incidence  angles  of  -10.9°  to  2.7°).  The  general 
layout  of  the  cascade  is  shown  in  Figure  27  on  page  52.  A  computer  code  that  generates 
the  coordinates  of  the  blade  and  a  summary  of  some  experimental  results  are  given  in 
Appendix  B. 

The  code  was  run  with  the  intermittency  constant  Gv  =  10.  Higher  values  of  Gy 
(above  100)  caused  numerical  problems  in  the  code.  The  onset  of  transition  was  first 
taken  at  the  point  where  it  was  observed  in  the  experiment.  At  the  lower  inlet  angle  it 
seems  that  a  better  agreement  with  the  experimental  results  can  be  obtained  by  delaying 
the  onset  of  transition  but  trying  to  implement  it  resulted  in  numerical  breakdown  of  the 
computation.  At  the  higher  inlet  angles,  better  agreement  with  the  experimental  results 
was  achieved  by  initiating  the  transition  earlier  (at  26%  chord  for  P  =  45.6°  and  at  21% 
for  P  =47.7°  as  compared  to  44%  and  3 6 To  chord  as  observed  in  the  experiment). 

1.  Displacement  Thickness 

Comparisons  of  the  experimental  data  to  the  computed  displacement  thickness 
are  shown  in  Figure  28  on  page  53  for  inlet  angle  of  34.1°,  in  Figure  29  on  page  54  for 
inlet  angle  of  36.3°,  in  Figure  30  on  page  55  for  inlet  angle  of  45.6°  and  in  Figure  31 
on  page  56  for  inlet  angle  of  47.7°. 

As  can  be  seen  in  the  figures,  there  is  a  good  agreement  between  the  actual  and 
the  computed  results  at  the  two  lower  angles  (/?  =  34.1°  and  /?  =  36.3°,  in  which  the 
incidence  angles  were  negative).  At  the  tw  o  higher  angles,  p  =  45.6°  and  />  =  47.7°  the 
computed  results  agree  with  the  actual  results  up  to  about  70%  chord,  and  then  the 
displacement  thickness  predicted  by  the  code  becomes  much  thicker  than  the  actual  one. 

The  code  predicted  a  large  flow  separation  area  starting  at  about  70%  chord  at 
the  lower  inlet  angles,  and  at  about  46%  chord  at  the  higher  inlet  angles.  This  flow 
separation  was  not  obsc-ved  in  the  experiment.  The  discrepancies  between  the  computed 
and  the  actual  results  behind  60%  to  70%  chord  can  be  explained  by  the  inaccurate 
calculations  by  the  code  due  to  the  large  separated  areas.  When  the  code  encounters 
separation,  several  approximations  are  made  (like  the  FLARE  approximation)  based  on 
the  assumption  that  the  separated  area  is  small.  When  the  separated  area  is  large,  these 
approximations  may  result  in  inaccurate  prediction  of  the  flow  field. 
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Figure  29.  C4  cascade  at  p  =  36.3°:  Displacement  thickness  comparison  with 
computed  results. 
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Fi  gure  30.  C4  cascade  at  (l  =  45.6°:  Displacement  thickness  comparison  with 
computed  results  (G=  10). 
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Figure  31.  C4  cascade  at  fl  =  47.7°:  Displacement  thickness  comparison  with 
computed  results  (G=  1U). 


2.  External  Velocity  and  Velocity  Profiles  Comparisons 

A  comparison  of  the  external  velocity  on  the  upper  surface  of  the  blade  is  shown 
in  figure  32  on  page  5S  for  inlet  angle  of  45.6°  and  in  Figure  33  on  page  59  for  inlet 
angle  of  4“.7\  It  can  be  seen  that  there  is  a  good  agreement  between  the  experimental 
and  the  computed  results  up  to  about  80"  c.  chord.  Near  the  trailing  edge  the  computed 
results  deviate  from  the  experimental  results  due  to  the  inaccuracy  in  the  calculations 
of  the  displacement  thickness. 

A  comparison  of  the  velocity  profiles  in  the  boundary  layer  at  50%  chord  is 
shown  in  Figure  34  on  page  60  for  inlet  angle  of  34.1°  and  in  Figure  35  on  page  61  for 
inlet  angle  of  36.3°.  The  agreement  between  the  calculated  velocity  profiles  and  the 
measured  velocity  profiles  is  very  good. 
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Figure  32.  C4  cascade  at  /J  =  45.6°:  External  velocity  distribution. 
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Figure  33.  C4  cascade  at  /?  =  47.7°:  External  velocity  distribution. 
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Figure  35.  C4  cascade  at  (1  =  36.3°:  Velocity  profile  at  50%  chord. 
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V.  CONCLUSIONS  AND  RECOMMENDATIONS 


A.  CONCLUSIONS 

The  interactive  viscous  inviscid  computer  code,  has  been  investigated  by  comparing 
its  predictions  of  boundary  layer  parameters  to  experimental  data. 

It  has  been  found  that  the  code  yields  reasonable  results  for  lightly  loaded  cascades, 
but  the  prediction  of  the  boundary  layer  thickness  on  the  suction  surface  of  highly 
loaded  cascades  deviates  significantly  from  experimentally  measured  data.  In  two  cases 
involving  highly  cambered  cascade  blades,  with  sharp  leading  edge,  the  code  failed  to 
run.  It  was  also  found  that  the  prediction  of  external  velocity  distribution  on  highly 
loaded  cascades  was  inaccurate. 

The  main  reasons  to  the  discrepancies  in  the  prediction  of  the  boundary  layer 
thickness  seem  to  be: 

1.  Inaccuracy  in  predicting  flow  parameters  in  regions  of  large  flow  separation  (due 
to  inadequate  transition  model  and  approximations  made  in  calculating  the  How  in 
separated  areas). 

2.  Inaccurate  turbulence  modelling. 

3.  Possible  violation  of  the  basic  assumptions  of  the  boundary  layer  theory  in  areas 
of  very  thick  boundary  layers. 

4.  The  wake  is  not  calculated  by  the  code.  The  result  is  inaccurate  flow  prediction 
near  the  trailing  edge. 

The  inaccuracy  in  the  prediction  of  the  external  velocity  distribution  in  highly  loaded 
cascades  is  due  to  the  interaction  law,  which  does  not  account  for  the  presence  of  adja¬ 
cent  blades. 


B.  RECOMMENDATIONS 

The  recommended  steps  in  order  to  improve  the  code  are: 

1.  Improving  the  interaction  law  by  assuming  a  distribution  of  sources  on  the  actual 
surface  (instead  of  the  assumption  of  a  Oat  plate),  letting  the  correction  term  to  the 
external  velocity  vary  across  the  boundary  layer  and  distributing  sources  on  the 
adjacent  blade  as  well  for  better  modelling  of  the  boundary  layer  effect  on  the  ex¬ 
ternal  velocity. 
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Changes  to  the  derivation  of  the  boundary  layer  equations  should  be  investigated 
to  allow  a  better  treatment  of  thick  boundary  layers  tlike  omitting  the  assumption 
of  ''/■  cy  =  0  across  the  boundary  layer). 

Different  turbulence  models  should  be  investigated. 

The  wake  should  be  included  in  the  calculations. 
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APPENDIX  A.  COMPUTER  CODE  LISTING 


>***************^****************************************************int00010 

VISCOUS-INVISCID  INTERACTION  PROGRAM  FOR  CASCADE  FLOWS  *****INT00020 

INT00040 

VERSION  3. A  INT00050 

JANUARY  87  INTO 0060 

INT00070 

THIS  VISCOUS-INVISCID  INTERACTION  METHOD,  CAPABLE  OF  COMPUTING  BOTH  INT00080 
SINGLE  AIRFOIL  AND  CASCADE  FLOWS,  WAS  DEVELOPED  BY  CEBECI  AND  INT00090 

COLLABORATEURS  AT  LONG  BEACH  STATE  AND  DOUGLAS  AIRCRAFT  COMPANY.  INT00100 
THE  CODE  APPLIES  TO  INCOMPRESSIBLE,  2-DIMENSIONAL,  STEADY  FLOWS  INT00110 
PAST  LINEAR,  ARBITRARILY  STAGGERED  CASCADES.  THE  METHODS  BASIC  INT00120 

INGREDIENTS  INC UDE  INT00130 

1.  A  FIRST  ORDL I  PANEL  METHOD  TO  SOLVE  LAPLACE'S  EQUATION,  INT00140 

2.  A  FINITE  DIFFERENCE  SCHEME  TO  SOLVE  THE  BOUNDARY  LAYER  EQUATIONS INTOO 150 


SUBJECT  TO  DIRECT  OR  INTERACTIVE  BOUNDARY  CONDITIONS,  INT00160 

3.  A  STRONG  INTERACTION  MODEL  TO  COUPLE  VISCOUS  AND  INVISCID  FLOW  INT00170 

RESULTS,  AND  INT00180 

4.  A  ZERO  EQUATION,  ALGEBRAIC  TURBULENCE  MODEL  TO  ESTIMATE  INTOO 190 

TURBULENT  SHEAR  STRESSES.  INT00200 

INT00210 

IN  SUMMARY,  THE  CODE  WILL  PROVIDE,  FOR  ATTACHED  AS  WELL  AS  MODERATE- I NT00220 
LY  SEPARATED  FLOWS  PAST  SINGLE  AIRFOILS  OR  CASCADES,  THE  FOLLOWING  INT00230 

1.  INVISCID  AND  VISCOUS  PRESSURE  DISTRIBUTIONS,  INT00240 


DISTRIBUTIONS  OF 

A.  LOCAL  SKIN  FRICTION  COEFFICIENT, 

B.  DISPLACEMENT  AND  MOMENTUM  THICKNESS,  AND 
VELOCITY  PROFILES  ACROSS  THE  BOUNDARY  LAYER. 


INT00230 

INT00240 

INT00250 

INT00260 

INT00270 

INT00280 

INT00290 

INT00300 

INT00310 


MODIFICATIONS  SINCE  VERSION  3.0:  INT00300 

1.  PRECISE  ASSIGNMENT  OF  BEGIN  OF  TRANSITION.  INT00310 

2.  CORRECTION  OF  AN  ERROR  IN  THE  CALCULATION  OF  MOMENTUM  THICKNESS.  INT00320 

3.  ADDITIONAL  PRINT  OPTION:  IP=-2  WILL  PROVIDE  AN  INPUT  FILE  (UNIT  INT00330 


NUMBER  12)  FOR  THE  PLOTTING  ROUTINE. 


COMMON  /BLCO/  NX , NXT , NP , NPT , NTR , IT , INVRS , NS , IP 
C  OMMON / B LOW / VN ( 100) 

COMMON/BLIN/  TITLE! 20) ,XC( 100) ,YC( 100) , ISG( 100) ,DELS( 100) , 

+  XCTR , XTR , ISTRP , I CYCLE , ICYTL ,XCTRS( 2) ,TRFIND( 2) 

COMMON/CASCDE/ INLET , SP .SINGLE , ALPHAA , ALPHAI , STAG 
COMMON/TRN/  PGAMTR , OMEGA , RTHETB , RTRANB 
COMMON/ PLOT/NVP( 2) ,NXVP( 20 , 2) , ICC 

DIMENSION  X0( 100) ,Y0( 100) ,X( 100) ,Y( 100) ,VCOM( 100) ,DLS( 100) , 
+  XS( 100) , YS( 100) ,XSTGR( 100) ,YSTGR( 100) ,DBPP(  100) 

DIMENSION  CASEID( 20) ,XCTRI( 2) , ITRI( 2) ,NBL(2) 

LOGICAL  SINGLE, TRFIND 
TRFIND( 1 )=  .FALSE. 

TRFIND( 2)=  .FALSE. 


INTO 0340 
INT00350 
INT00360 
INT00370 
INT00380 
INT00390 
INT00400 
INT00410 
INT00420 
INT00430 
INT00440 
INT00450 
INT00460 
INT00470 
INT00480 
INT00490 
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ICASE  =  0 

READ(  5 ,5 ,END=999)  TITLE 
FORMAT! 20A4) 

ICASE  =  ICASE  +  1 
REWIND  3 
READ  (5,10) 

FORMAT! IX) 

READ  (5,20)  ITRI(l) ,ITRI(2) , IRST, ICYTL, IP 
FCRMAT( 1615 ) 

READ  (5,10) 

READ  (5,25)  INLET, ISTAG.ALPHAI .STAG, SP.PGAMTR, OMEGA 
FORMAT! 215 , 5F10.  0) 

READ  ( 5 , 2 7 )RN , XCTRI ( 1 ) , XCTRI ( 2 ) , ALPHAA 
FORMAT! 4E 10.  0) 

IF  (IP. EQ. -2)  THEN 

READ  (5,20)  NVP( 1) ,NVP( 2) 

IF  (NVP(l).NE.O)  READ  (5,20)  (NXVP( I , 1) , 1=1 ,NVP( 1) ) 
IF  (NVP( 2) . NE. 0)  READ  (5,20)  (NXVP( I ,2) , 1=1 ,NVP(2) ) 
END  IF 

IF  (ICASE  .EQ.  1)  READ  (5,20)  N,NI 
I READ  =  1 
I BLOW  =  1 
SINGLE  =  .FALSE. 

IF  (SP  .LE.  0.0)  SINGLE  =  .TRUE. 

N  =  N  -  1 
Nl=  N  +  1 

IF  (ICASE  .GT.  1)  THEN 
N1  =  N1SAVE 
N  =  N  -  1 
GOTO  53 
END  IF 


WRITE( 12,20)  N+l ,NVP( 1 ) ,NVP( 2) ,90,70, INLET 
IF  (NVP(l).NE. 0)  WRITE( 12,20)  (NXVP( I , 1) , 1=1 ,NVP( 1) ) 
IF  ( NVP( 2) . NE. 0)  WRITE( 12,20)  (NXVP( I , 2) , 1=1 ,NVP( 2) ) 
IF  ( INLET.  NE.  1)  WRITE( 12,80)  RN, ALPHAA 
IF  ( INLET. EQ. 1)  WRITE(12,80)  RN, ALPHAI 
WRITE! 12,82)  (X0( I ) , 1=1 , N+l) 

WRITE( 12,82)  (Y0( I) , 1=1, N+l) 

FORMAT! 2E 15.  5) 

F0RMAT( 8F10.  6) 

END  IF 


INT00500 
INT005 10 
INT00520 
INT00530 
INT00540 
INT00550 
INT00560 
INT00570 
INT00580 
INT00590 
INT00600 
INT00610 
INT00620 
INT00630 
INT00640 
INT00650 
INT00660 
INT00670 
INTO  0680- 
INT00690 
INT00700 
INT00710 
INT00720 
INT00730 
INT00740 
INT00750 
INT00760 
INT00770 
INT00780 
INT00790 
INT00800 
INT00810 


c 

READ  (5,10) 

INT00820 

READ  (5 , 3C )  (X0( I) , 

Y0( I ) ,  1=1, N+l) 

INT00830 

30 

F0RMAT( 2F10.  0) 

INT00840 

GO  TO  50 

INT00850 

C 

INT00S60 

C40 

READ(5,10) 

INT00870 

40 

READ( 5,45)  (XO(I)  , 

1=1, N+l) 

INT00880 

C 

READ! 5,10) 

INT00890 

READ(5 ,45)  (YO(I)  , 

1=1, N+l) 

INT00900 

45 

FORMAT ( 6F10.  0) 

INT00910 

C 

INT00920 

50 

CONTINUE 

INT00930 

IF  ( IP. EQ. -2)  THEN 

INT00940 

INT00950 
INT00960 
INT00970 
INT00980 
INT00990 
INT01000 
INT01010 
INT01020 
INTO  1030 
INTO  1040 
INTO  1050 
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non 


NRITE  =  (NI+l)/2 

IMIS  =  (Sl-l)/2+l 

IF(  (N'l/2*2)  .EQ.  Nl)  IMIN  =  Nl/2 

CALL  TRGRID  (Nl ,XO,YO,NI .NRITE, 0. 5 , IMIN, RAD, 1.NXSS1) 

N1SAVE  =  Nl 
S3  rONTTNCF 

ALPHAA  =  0.0174533  *  ALPHAA 
ALPHAI  =  0.0174533  *  ALPHAI 
STAG  =  0.0174533  *  STAG 
C 

IF  (INLET  .EQ.  0)  THEN 

ALPHA  =  ALPHAA 

ELSE 

ALPHA  =  ALPHAI 
END  IF 
C 

IF  (ISTAG  .NE.  0)  THEN 

CALL  STAGRC N 1 , STAG , XO , YO , XSTGR , YSTGR) 

ELSE 

DO  55  I  =  1  ,  Nl 
XSTGR(I)  =  XO(I) 

YSTGR(I)  =  YO( I ) 

55  CONTINUE 

END  IF 

READ  DATA  FROM  VISCOUS  CAL. 

I CYCLE  =  0 
60  ICYCLE  =  ICYCLE  +  1 

C 

CALL  POTNL( N 1 , IRST , ALPHA , CHORD , XO , YO , XSTGR , YSTGR , X , Y , DLS , VCOM , 

+  DBPP) 

IF  (ICYCLE  .GT.  ICYTL)  THEN 
REWIND  3 

WRITE  ( 3)N1 , (XO( I ) , YO( I ) ,DLS( I) , VN( I) ,DBPP( I) , 1=1 ,N1) 

GOTO  1 
END  IF 
C 

IF  (ISTAG  .NE.  0)  THEN 

DO  70  I  =  1  ,  Nl- 1 

X(I)  =  0.5  *  (XO(I)+XO(I+l)) 

Y( I )  =0.5  *  (YO( I)+YO( 1+1) ) 

70  CONTINUE 
END  IF 
C 

CALL  CASB LP( N 1 , XO , YO , X , Y , XS , YS , DLS , VCOM , DBPP , RN 
+  ,NBL,ITRI,XCTRI, TITLE) 

C-0  TO  60 
999  CONTINUE 
STOP 
END 
C 

SUBROUTINE  POTNL( Nl , IRST , ALPHA , CHORD , XO , YO , XSTGR , YSTGR , X , Y , DLS , 
+  VCOM, DBPP) 

C 

COMMON/BLOW/VN( 100) 


INTO  1060 
INT01070 
INT01080 
INTO  1090 
INTO 1100 
INT01110 
INTO 1120 
INT01130 
INTO  1140 
INTO  1150 
INT01160 
INT01170 
INTO  1180 
INT01190 
INT01200 
INT01210 
INT01220 
INTO  1230 
INT01240 
INT01250 
INTO 1260 
INT01270 
INTO 1280 
INTO 1290 
INT01300 
INTO  13 10 
INT01320 
INTO 1330 
INTO 1340 
INTO 1350 
INTO 1360 
INTO 1370 
INTO 1380 
INT01390 
INTO 1400 
INTO 1410 
INTO 1420 
INTO 1430 
INTO  1440 
INTO 1450 
INT01460 
INTO  1470 
INT01480 
INT01490 
INT01500 
INT015 10 
INTO  1520 
INT01530 
INT01540 
INT01550 
INTO 15 60 
INT01570 
INTO 15 80 
INTO  1590 
INT01600 
INT01610 
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COMMON  /BLCO/  NX ,NXT,NP ,NPT,NTR , IT , INVRS ,NS , IP 
C0MM0N/BLIN/TITLE(20) ,XC( 100) ,YC( 100) , ISG( 100) ,DELSC 100) , 

+  XCTR , XTR , ISYRP , I CYCLE , ICYTL, XCTRS (2) ,TRFIND( 2) 

COMMON/ CASCDE  -  INLET , SP , S INGLE , ALPHAA , ALPHAI , STAG 
C  SIMPLE  SOURCE  POTENTIAL  CODE 

DIMENSION  AQFF( 100 , 100) ,  EOFF( 100 , 100) ,  XP(100),  YP( 100) ,X( 100) , 
+  S(100),Cf 100),  D( 100) ,VTAN(3 , 100) , VNOR( 3, 100), R( 3, 100), 

+  VC CM( 100) ,SIGCOM( 100) ,CP( 100) ,XO(100) ,Y0(  100) 

+  , VNC( 100) , D 1 C 100) ,D2( 100) ,D3( 100) ,S0( 100) ,SC( 100) 

+  ,XOFF( 100) , YOFF( 100) ,T( 3 , 100) ,VTCOM( 100) ,VNCOM( 100) 

+  ,XS( 100) ,YS(  100) ,SOFF( 100) ,COFF( 100) ,XPOFF( 100) ,YPOFF( 100) 

+  ,Y(  100) ,SIG(3, 100) ,DLS( 100) ,DLSC( 100) ,A( 100,100) ,B(  100,100) 

+  , XSTGR( 100) ,YSTGR( 100) ,VUT(3) ,VLT(3) ,VUN( 3) ,VLN( 3) ,DBPP( 100) 

+  ,CPI(  100) ,XOS( 100) ,YOS( 100) ,DBPPC( 100) 

REAL  NUM1  ,  NUM2 
LOGICAL  OFF, SINGLE 
OFF  =  .FALSE. 

PI  =  3. 1415S2 
CM  =  0.  0 
N  =  N1  -  1 

IF  (I CYCLE  .EQ.  1)  THEN 
IF  (IRST  .EQ.  0)  THEN 
DO  10  1=1, N1 
CLS(I)  =  0.  0 
VN  (I)  =  0.0 
DBPP( I )=  0.  0 
10  CONTINUE 

ELSE 

DO  5  I  =  1  ,  N1 

XOS(I)  =  XO(I) 

YOS(I)  =  YO(I) 

5  CONTINUE 

READ  (3)  NT, (XS( I ) , YS( I) ,DLSC( I) , VNC( I ) ,DBPPC( I) , 1=1 , NT) 

XMIN  =  XS(1) 

DO  15  I  =  2  ,  NT 

IF  (XS(I)  . GT.  XMIN)  GOTO  15 

XMIN  =  XS(I) 

IMIN  =  I 

15  CONTINUE 

DO  17  I  =  1  ,  NT 
IF  (I  .  LT.  IMIN)  GOTO  16 
XS(I)  =  XS(I)  -  XMIN 
GOTO  17 

16  XS(I)  =  XMIN  -  XS(I) 

17  CONTINUE 
C 

XMIN  =  XOSCI) 

DO  20  I  =  2  ,  N1 

IF  (XOS(I)  .GT.  XMIN)  GOTO  20 

XMIN  =  XOS(I) 

IMIN  =  I 
20  CONTINUE 

DO  22  I  =  1  ,  N1 

IF  (I  . LT.  IMIN)  GOTO  21 

XOS(I)  =  XOS(I)  -  XMIN 


INT01620 
INTO  1630 
INTO  1640 
INT01650 
INT01660 
INT01670 
INT01680 
INTO 1690 
INTO  1700 
INTO  17 10 
INT01720 
INTO  1730 
INT01740 
INTO  1750 
INTO 1760 
INT01770 
INTO  17 80 
INTO  1790 
INT01800 
INTO  18 10 
INTO  1820 
INT01830 
INT01840 
INTO  1850 
INTO  1860 
INTO  1870 
INT01880 
INTO  1890 
INTO  1900 
INT01910 
INT01920 
INT01930 
INT01940 
INTO  1950 
INT01960 
INTO 19 70 
INT01980 
INT01990 
INT02000 
1NT02010 
INT02020 
INT02030 
INT02040 
INT02050 
INT02060 
INT02070 
INT02080 
INT02090 
INT02100 
INT02110 
INT02120 
INT02130 
INTO 2 140 
INT02150 
INT02160 
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GOTO  22 

INT02170 

21 

XOS(I)  =  XMIN  -  XOS(I) 

INT02180 

22 

CONTINUE 

INT02190 

C 

INT02200 

CALL  DIFF3(NT,XS ,DLSC ,D1 ,D2 ,D3 ,0) 

INT02210 

CALL  INTRP3(NT,XS ,DLSC ,D1 ,D2 ,D3 ,N1 ,XOS ,DLS) 

INT02220 

CALL  AMEAN  ( 1 ,N1 ,XOS ,DLS , 1) 

INT02230 

CALL  DIFF3  (NT,XS , VNC ,D1 ,D2 ,D3 , 0) 

INT02240 

CALL  I NTRP3 ( NT , XS , VNC , D 1 , D2 , D3 , N 1 , XOS , VN ) 

INT02250 

CALL  AMEAN  (1,N1,X0S,VN,1) 

INT02260 

CALL  DIFF3  (NT,XS}DBPPC,D1,D2,D3,0) 

INT02270 

CALL  INTRP3( NT , XS , DBPPC , D 1 , D2 , D3 , N1 , XOS , DBPP ) 

INT02280 

CALL  AMEAN  ( 1, Nl, XOS, DBPP, 1) 

INT02290 

END  IF 

INT02300 

END  IF 

INT02310 

DO  30  1=1, Nl 

INT02320 

XP(I)  =  XSTGR(I) 

INT02330 

YP( I)  =  YSTGR(T) 

INT02340 

30 

CONTINUE 

INT02350 

C 

CALCULATE  GEOMETRIC  QUANTITIES 

INT02360 

DO  100  J=1  ,N 

INT02370 

VNC(J)  =  0.5  *  (VN( J)  +  VN( J+l) ) 

INT02380 

X(J)=  . 5*(XP(J)+XP(J+1)) 

INT02390 

Y(J)=  .  5*(YP(J)+YP(J+1)) 

INT02400 

D(  J)  =  SQRT((XP(J+1)-XP(J»**2  +  (  YP(  J+l)  -YP(  J)  )**2) 

INTO  24 10 

C(J)=  (XP( J+l) -XP( J) )/D( J) 

INT02420 

S(J)=  (YP(J+1)-YP(J))/D(J) 

INT02430 

100 

CONTINUE 

INT02440 

C 

INT02450 

IF  (  INLET  .NE.  0  .AND.  .NOT.  SINGLE)  THEN 

INT02460 

SUM  =  D( 1) 

INT02470 

DO  35  J  =  2,  N 

INT02480 

SUM  =  SUM  +  D( J) 

INT02490 

35 

CONTINUE 

INT02500 

Q  =  2.  0  *  PI  *  SUM  /  SP 

INT02510 

ELSE 

INT02520 

Q  =  0.  0 

INT02530 

END  IF 

INT02540 

C 

INT02550 

C 

CALCULATE  NORMAL  AND  TANGENTIAL  MATRICES 

INT02560 

102 

CONTINUE 

INT02570 

IF  (SINGLE)  THEN 

INT02580 

IF  (  .NOT.  OFF)  THEN 

INT02590 

DO  120  1=1, N 

INT02600 

DO  110  J=1 ,N 

INT02610 

IF  (J  . EQ.  I)  GO  TO  105 

INT02620 

XX=  (X( I ) -X( J) )*C( J)  +  (Y( I) -Y( J) )*S( J) 

INT02630 

YY=-(X( I ) -X( J) )*S( J)  +  (Y(I)-Y(J) )*C( J) 

INT02640 

UU=  LOG( ( (XX+.  5*D ( J ) ) ** 2+YY** 2 ) / ( ( XX - .  5*D( J) )**2+YY**2) ) 

INT02650 

VV=  2. *ATAN2(YY*D(J),  XX** 2+YY** 2 - ( . 5*D( J ) )**2 ) 

INT02660 

SS=  S(I)*C(J)  -  C(I)*S(J) 

INT02670 

CC=  C(I)*C(J)  +  S(I)*S(J) 

INT02680 

A( I , J)=  -UU*SS  +  VV*CC 

INT02690 

B(I,J)=  UU*CC  +  VV*SS 

INT02700 

GO  TO  110 

INT027 10 

105 

A( I , J)  =  6. 2831853 

INT02720 
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B( I , J)  =  0. 0 
110  CONTINUE 
120  CONTINUE 
ELSE 

DO  149  1  =  1,  N 
DO  130  J=1,N 

XX=  (XOFF(I)-X(J))*C(J)  +  ( YOFF( I) -Y( J) )*S( J) 

YY=-( XOFFC I) -X( J) )*S( J)  +  (YOFF( I) -Y( J) )*C( J) 

UU=  LOG( ( (XX+.  5*D( J) )**2+YY**2)/( (XX-.  S*D( J) )**2+YY-  *2)) 
VV=  2. * AT AN 2 ( Y Y*D ( J ) ,  XX**2+YY**2-( .  5*D( J) )**2) 

SS=  SOFF( I)*C( J)  -  COFF( I)*S( J) 

CC=  COFF( I )*C( J)  +  SOFF( I )*S( J) 

AOFF(  I ,  J)=  -UU*SS  +  W*CC 
BOFF(I,J)=  UU*CC  +  VV*SS 
130  CONTINUE 
140  CONTINUE 
END  IF 
ELSE 

IF  (  .NOT.  OFF)  THEN 

DO  50  1=1, N 

DO  40  J=1 ,N 

IF  (J  .EQ.  I)  GO  TO  45 

XX=  (X(I)-X(J))*C(J)  +  ( Y( I ) -Y( J) )*S( J) 

YY=-(X( I) -X( J) )*S( J)  +  (Y(I)-Y(J) )*C( J) 

DX1  =  PI  *(X(I)-XP(J))  /  SP 

DY1  =  PI  *(Y( : ) -YP( J) )  /  SP 

DX2  =  PI  *(X(I)-XP(J+1))  /  SP 

DY2  =  PI  *( Y( I ) -YP( J+l) )  /  SP 

R1SQ  =  (C0SHfTv:i))**2  -  (COS(DYl) )**2 

R2S0  =  (C0SH(DX2))**2  -  ( C0S(DY2) )**2 

UU  =  L0G(R1SQ/R2SQ) 

NUM1  =  DX1  *  COSH(DXl)  *  SIN(DYl)  - 

+  DY1  *  SINH(DXl)  *  COS(DYl) 

DNUM1=  DX1  *  SINH(DXl)  *  COS(DYl)  + 

+  DY1  *  COSH(DXl)  *  SIN(DYl) 

NUM2  =  DX2  *  COSH(DX2)  *  SIN(DY2)  - 

+  DY2  *  SINH(DX2)  *  COS(DY2) 

DNUM2=  DX2  *  SINH(DX2)  *  COS(DY2)  + 

+  DY2  *  COSHCDX2)  *  SIN(DY2) 

EXV  =  2.0  *  ATAN2 ( NUM2 , DNUM2 )  -  2. 0  *  ATAN2 ( NUM 1 , DNUM 1 ) 
VV=  2. *ATAN2(YY*D(J) ,  XX**2+YY**2-(.  5*D(J))**2) 

VV  =  VV  +  EXV 

SS=  S(I)*C(J)  -  C(I)*S(J) 

CC=  C(I)*C(J)  +  S(I)*S(J) 

A(I,J)=  -UU*SS  +  VV*CC 
B( I , J)=  UU*CC  +  VV*SS 
GO  TO  40 

45  A( I , J)  =  6. 2831853 

B(  I ,  J)  =  0.0 

40  CONTINUE 

50  CONTINUE 

ELSE 

DO  70  1=1, N 
DO  60  J=1 ,N 

XX=  (XOFF( I ) *X( J) )*C( J)  +  (YOFF(I)-Y(J))*S(J) 

YY=-(XGFF( I ) -X( J) )*S( J)  +  ( YOFF( I ) -Y( J) )*C( J) 


INT02730 
INT02740 
INT02750 
INT02760 
INT02770 
INT02780 
INT02790 
INT02800 
INT028 10 
INT02820 
INT02830 
INT02840 
INT02850 
INT02860 
INT02870 
INT02880 
INT02890 
INT02900 
INT029 10 
INT02920 
INT02930 
INT02940 
INT02950 
INT02960 
INT02970 
INT02980 
INT02990 
INT03000 
INT03010 
INT03020 
INT03030 
INT03040 
INT03050 
INT03060 
INT03070 
INT03080 
INT03090 
INT03100 
INT03110 
INT03120 
INT03130 
INT03140 
INT03150 
INT03160 
INT03170 
INT03180 
INT03190 
INT03200 
INT03210 
INT03220 
INT03230 
INT03240 
INT03250 
INT03260 
INT03270 
INT03280 
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DX1  =  PI  *(XOFFf I)-XP(J))  /  SP 

INT03290 

DY1  =  PI  *( YOFF( I ) -YP( J) )  /  SP 

INTO 3 3 00 

DX2  =  PI  *(X0FF(I)-XP(J+I))  /  SP 

INT03310 

DY2  =  PI  *(  Y'0FF(  I )  -YP(  J+l ) )  /  SP 

INT03320 

R1SQ  =  ( COSHC DX 1 ) ) ** 2  -  (COS(DYl) )**2 

INT03330 

R2SQ  =  (COSH(DX2))**2  -  ( C0S(DY2) )**2 

INT03340 

UU  =  L0G(R1SQ/R2SQ) 

INT03350 

NUM1  =  DX1  *  COSH(DXl)  *  SIN(DYl)  - 

INT03360 

+  DY1  *  SINK(DXl)  *  COS(DYl) 

INT03370 

DNUM1=  DX1  *  SINH(DXl)  *  COS(DYl)  + 

INT03380 

+  DY1  *  COSH(DXl)  *  SIN(DYl) 

INT03390 

NUM2  =  DX2  *  COSH( DX2)  *  SIN(DY2)  - 

INT03400 

+  DY2  *  SINHCDX2)  *  COS(DY2) 

INT03410 

DNUM2=  DX2  *  SINH(DX2)  *  C0S(DY2)  + 

INT03420 

+  DY2  *  COSH( DX2)  *  SINCDY2) 

INT03430 

EXV  =  2.0  *  ATAN2C NUM2 , DNUM2 )  -  2.  0  *  AT AN 2 ( NUM 1 , DNUM 1 ) 

INT03440 

VV=  2.  *ATAN2(YY*D( J) ,  XX**2+YY**2-( • 5*D( J) )**2) 

INT03450 

VV  =  VV  +  EXV 

INT03460 

SS=  SOFFC I )*C( J)  -  COFFCI)*SCJ) 

INT03470 

CC=  C0FF(I)*C(J)  +  SOFFC  I )*S( J) 

INT034S0 

AOFF(I,J)=  -UU*SS  +  VV*CC 

INT03490 

BOFF(  I ,  J)=  UU*CC  +  W*SS 

INT03500 

60 

CONTINUE 

INT03510 

70 

CONTINUE 

INT03520 

END  IF 

INT03530 

END  IF 

INT03540 

C 

NORMAL  AND  TANGENTIAL  COMPONENTS  OF  FUNDAMENTAL  SOLUTIONS 

INT03550 

c 

INT03560 

DO  160  1=1, N 

INT03570 

SUMR=  0. 

INT03580 

SUMT=  0. 

INTO 35 90 

IF  (  .NOT.  OFF)  THEN 

INT03600 

RC1,I)=  S( I )+VNC( I ) /COS (ALPHA) 

INT03610 

TC 1 , 1)=  CCD 

INT03620 

R(  2 , 1 )=  -C(I) 

INT03630 

T(  2 , 1 )=  S(I) 

INT03640 

DO  145  J=1 ,N 

INT03650 

SUMR  =  SUMR  +  B( I , J) 

INT03660 

SUMT  =  SUMT  +  ACI.J) 

INT03670 

145 

CONTINUE 

INT03680 

ELSE 

INT03690 

R( 1, I)  =  SOFFC I) 

INT03700 

TC 1,1)  =  COFF(I) 

INT037 10 

R( 2 , I )  =-COFF( I ) 

INT03720 

T( 2 , I )  =  SOFFC I) 

INT03730 

DO  150  J=1,N 

INT03740 

SUMR=  SUMR  +  BOFF(I.J) 

INT03750 

SUMT=  SUMT  +  A0FF( I , J) 

INT03760 

150 

CONTINUE 

INT03770 

END  IF 

INT03780 

R( 3 , I )=  SUMR 

INT03790 

T( 3 , I )=  SUMT 

INTO 3800 

160 

CONTINUE 

INT03810 

C 

INT03S20 

IF  (  OFF  )  GO  TO  275 

INT03830 

C 

DECOMPOSITION  OF  MATRIX  A 

1NT03840 

70 


DO  230  1=1, S-l 

INT03850 

DO  220  K=I+1 , N 

INT03860 

A(K, I )=  A(K, I)/A( I , I) 

INT03870 

DO  210  J=I+1 ,N 

INT03880 

A(K, J)=  A(K, J) -  A(K, I)*A( I , J) 

INT03890 

210 

CONTINUE 

INT03900 

220 

CONTINUE 

INT039 10 

230 

CONTINUE 

INT03920 

C  OPERATE  ON  FUNDAMENTAL  RIGHT  SIDES  WITH  LOWER  TRIANGULAR 

INT03930 

DO  270  K=1 , 3 

INT03940 

DO  260  J=1,N-1 

INT03950 

DO  250  I=J+1 ,N 

INT03960 

R(K, I )=  R(K, I)  -  A( I , J)*R(K, J) 

INT03970 

250 

CONTINUE 

INT03980 

260 

CONTINUE 

INT03990 

270 

CONTINUE 

INT04Q00 

C  ] 

BACK  SOLUTION 

INT04010 

DO  300  K=1 , 3 

INT04020 

DO  290  I =N , 1 , - 1 

INT04030 

SUM=  0. 

INT04040 

DO  280  J=N , 1+1 , - 1 

INT04050 

SUM=  SUM  +  A( I , J)*SIG(K, J) 

INT04060 

2S0 

CONTINUE 

INT04070 

SIG( K , I )=  (R(K, I ) -SUM) / A( 1,1) 

INT04080 

290 

CONTINUE 

INT04090 

300 

CONTINUE 

INT04100 

OFF  =  .TRUE. 

INT04110 

C 

INT04120 

c 

ADD  DIS-PLACE  VERTICALLY  TO  THE  BODY  TO  GENERATE 

INT04130 

c 

DISPLACEMENT  SURFACE 

INT04140 

c 

INT04150 

DO  305  I  =  2  ,  N 

INT04160 

SOFF(I)  =  (S(I)*D(I-1)+S(I-1)*D(I))/(D(I-1)+D(I)) 

INT04170 

COFFCI)  =  ( C( I)*D( I-l)+C( I-1)*D( I) )/(D( I -1)+D( I) ) 

INT04180 

305 

CONTINUE 

INT04190 

SOFF(l)  =  2. 0*S(1)  -  SOFF( 2) 

INT04200 

S0FF(N1)=  2. 0*S(N)  -  SOFF(N) 

INT04210 

COFF(l )  =  2. 0*C( 1)  -  COFF( 2) 

INT04220 

C0FF(N1)=  2. 0*C(N)  -  COFF(N) 

INT04230 

DO  306  I  =  1  ,  N1 

INT04240 

XPOFF(I)  =  XP(I)  -  SOFF(I)  *  DLS(I) 

1NT04250 

YPOFF(I)  =  YP( I )  +  COFF(I)  *  DLS(I) 

INT04260 

306 

CONTINUE 

INT04270 

DO  307  I  =  1  ,  N 

INT04280 

XOFF(I)  =  0.5  *  (XPOFF(I)  +  XPOFF(I+l)) 

INT04290 

YOFF(I)  =  0.5  *  (YPOFF(I)  +  YPOFF(I+l)) 

INT04300 

DOFF  =  SQRT( (TPOFFC 1+1 ) -XPOFF( I ) )**2  + 

INT04310 

+  QPOFF(  1+1 )  -YPOFF(  I )  )**2  ) 

INT04320 

COFF(I)  =  (XPOFF( 1+1) -XPOFFC I ) )/DOFF 

INT04330 

SOFF(I)  =  (YPO^-( I+l)-YPOFF(I))/DOFF 

INT04340 

307 

CONTINUE 

INT04350 

GO  TO  102 

INT04360 

C  ' 

CALCULATION  OF  SURFACE  VELOCITIES  FOR  THE  FUNDAMENTAL  SOLUTIONS 

INT04370 

C 

INT04380 

275 

DO  330  K=1 , 3 

INT04390 

DO  320  1=1, N 

INT04400 

71 


SUMT=  T(K, I) 

SUMN--RC K ,  I ) 

DO  310  j=l  ,N 

SUMT=  SUMT  +  BOFF( I , J)*SIG(K, J) 

SUMN=  SUMS  +  AOFF(I,J)*SIC-(K,J) 

310  CONTINUE 

VTAN ( K , I ) =  SUMT 
VNOR(K, I)=  SUMN 
320  CONTINUE 
C 

DOFF1  =  SQRT( (XPOFF( 2) -XPOFF( 1) )**2+( YPOFF( 2) -YPOFF( 1) )**2) 
DOFF 2  =  SQRT( (XPOFF( 3) -XPOFF( 2) )**2+( YPOFF( 3) -YPOFF( 2) )**2) 
DOFFN  =  SQRT( ( XPOFF( N+l ) -XPOFF(N) )**2+( YPOFF( N+l ) - YPOFF(N) ) 

+  **2) 

DOFFN 1=  SQRTC (XPOFF(N) -XPOFF(N-l) )**2+( YPOFF(N) -YPOFF(N-l) ) 

+  **2) 

VTJT(K)  =  VTAN(K,N)  +  DOFFN  *  ( VTAN(K,N) -VTAN(K,N-1) )/ 

+  ( DOFFN+DOFFN 1 ) 

VUN(K)  =  VNOR(K,N)  +  DOFFN  *  (VNOR(K,N)-VNOR(K,N-l))/ 

+  ( DOFFN+DOFFN 1) 

VLT(K)  =  VTAN(K, 1)  +  DOFF1  *  ( VTAN(K , 1) -VTAN(K, 2) )/ 

+  (DOFF1+DOFF2) 

VLN(K)  =  VNOR(K, 1)  +  DOFF1  *  ( VNOR(K, 1) -VNOR(K, 2) V 
+  (DOFF1+DOFF2) 

330  CONTINUE 

C  OUTPUT  FUNDAMENTAL  SOLUTIONS 

IF  (ICYCLE  .EQ.  1  .OR.  ICYCLE  .  GE.  ICYTL-1  .OR.  IP  .  GE.  0) 

+  WRITE (6,335)  TITLE 
335  F0RMAT( 1H1 , ///20A4//) 

C  DO  360  K=1 , 3 

C  WRITE  (6,340)  K 

C340  FORMAT( / /// , 1H  FUNDAMENTAL  SOLUTION  NUMBER  ',12////) 

C  WRITE( 6 , 345 ) 

345  FORMAT( 3X, ' I ' ,8X,'X' ,11X,'Y' ,10X,’VT’ ,10X,'VN' ,8X,'SIG'  ///) 
C  WRITE( 6 , 375 )  1,  XP(1)  ,  YP(1),  VLT(K) ,VLN(K) 

C  DO  350  1=1, N 

C  WRITE( 6 , 375 )  I  ,  X(I),  Y(I),  VTAN(K,I) ,VNOR(K,I) ,SIG(K,I) 

C350  CONTINUE 

C  WRITE( 6 , 375)  N1  ,  XP(N1) ,YP(N1) , VUT(K) ,VUN(K) 

C360  CONTINUE 

C  COMBINED  FLOW  AT  ANGLE  OF  ATTACK 
C 

IF  (  INLET  . NE.  0)  THEN 

YYY  =  ( ( VUT( 3)+VLT( 3) )*TAN( ALPHAI )+( VUT( 1)+VLT( 1 ) )*Q)/ 

+  ( ( VUT( 3 ) +VLT( 3) ) -( VUT( 2) +VLT( 2 ) )*Q ) 

XXX  =  -( ( VUT( 1)+VLT( 1) )+( VUT(2)+VLT(2) )*TAN( ALPHAI) )/ 

+  ( ( VUT( 3 ) +VLT( 3) ) -( VUT( 2 ) +VLT( 2 ) )*Q) 

ALPHA  =  ACOS( 1. 0/SQRT( 1.  0+YYY**2) ) 

COSAL  =  COS( ALPHA) 

SINAL  =  SIN( ALPHA) 

W  =  XXX/SQRT( 1.  0+YYY**2) 

ELSE 

COSAL  =  COS( ALPHA) 

SINAL  =  SIN( ALPHA) 

W=-((VLT(  1)+VUT(  1)  ),vCOSAL+(  VLT(  2)+VUT(  2)  ),VSINAL)  / 

+  ( VLT( 3)+VUT( 3) ) 


INT04410 
INT04420 
INTO 44 30 
INT04440 
INT04450 
INT04460 
INT04470 
INT04480 
INT04490 
INTU4500 
INT045 10 
INT04520 
INT04530 
INT04540 
INT04550 
INT04560 
INT04570 
INT04580 
INT04590 
INT04600 
INT04610 
INT04620 
INT04630 
INT04640 
INT04650 
INT04660 
INT04670 
INT04680 
INT04690 
INT04700 
INT047 10 
INT04720 
INT04730 
INT04740 
INT04750 
INT04760 
INT04770 
INT04780 
INT04790 
INT04800 
INT04810 
INT04820 
INT04830 
INT04840 
INT04850 
INT04860 
INT04870 
INT04880 
INT04890 
INT04900 
INT049 10 
INT04920 
INT04930 
INT04940 
INT04950 
INT04960 


o  n  n 


END  IF 

C  FORCE  COEFFICIENT  CALCULATION 
SUM1=  0. 

SUMX=  0. 

SUMY=  0. 

DO  390  1=1, N 
SUM1=  SUM1+  D( I ) 

SUMX=  SUMX-  VCOM( I )**2*S( I )*D( I ) 
SUMY=  SUMY+  VCOM( I )**2*C( I )*D( I ) 
390  CONTINUE 
C  FIND  MAN.  CHORD  LENGTH 
XOMIN  =  X0(1) 

DO  395  I  =  2  ,  N1 

IF  (  XO(I)  . GT.  XOMIN)  GOTO  395 

XOMIN  =  XO(I) 

395  CONTINUE 

CHORD  =  XO(Nl)  -  XOMIN 

CL1=  SUM 1*25 . 1327  4*W / CHORD 

CL2=  ( SUMY"*COSAL-SUMX*S  INAL) / CHORD 

CD  =  ( SUMX*CCSAL+SUMY*SINAL) /CHORD 
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INT05170 

CALCULATING  PARAMETERS  FOR  INLET  VELOCITY  AS  MODULUS  OF  NOMORIZED  VELINT05180 

INT05190 

IF  (.NOT.  SINGLE)  I TEN  INT05200 

NUM 1  =  SIN(ALPHA)+CLl*CH0RD/(4. 0*SP)  INT05210 

ALPHID  =  ATAN2 ( NUM 1, COS (ALPHA))  INT05220 

NUM 1  =  S IN( ALPHA) -CLl*CHORD/( 4. 0*SP)  INT05230 

ALPHED  =  AT AN2 ( NUM 1, COS (ALPHA))  INT05240 

NUM 1  =  CLl*CHORD/( 2.  0*SP)*COS( ALPHA)  INT05250 

DNUM1=  1. 0-(CLl*CH0RD/(4. 0*SP))**2  INT05260 

DALPKA  =  ATAN2 ( NUM 1 , DNUM 1 )  INT05270 

UOUI  =  (TAN( ALPHID) -TAN(ALPHED) )*( 2.  0*SP/CHORD*COS( ALPHID) ) /CL1  INT05280 
CLI  =  CL1  *  UOUI** 2  INT05290 

UIOU  =  1.  O/UOUI  INT05300 

VEXIT  =  COS ( ALPHA ) / COS ( ALPHED )  INT05310 

ELSE  INT05320 


ALPHID  =  ALPHA 
ALPHED  =  ALPHA 
DALPHA  =0.0 
UOUI  =  1.0 
UIOU  =1.0 
CLI  =  CLI 
VEXIT  =1.0 
END  IF 

FAC  =  180. 0/PI 
ALPHID  =  ALPHID  *  FAC 
ALPHED  =  ALPHED  *  FAC 
DALPHA  =  DALPHA  *  FAC 
ALPHAD  =  ALPHA  *  FAC 

IF  (ICYCLE  .EQ.  1  .OR.  ICYCLE  . GE.  ICYTL-1  .OR.  IP  .CE.  0)  THEN 

WRITE(6,370)  ALPHAD  , ALPHID .ALPHED, DALPHA, UIOU, VEXIT 

FORMAT  (////, 1H  ,  ’COMBINED  FLOW  AT  AVERAGE  ANGLE  OF  ATTACK  * 

+  F8.3,  4X, ’DEGREES’ ,  / , 1H  ,17X,’ INLET  ANGLE  OF  ’, 

+  ’ATTACK  =  ’ ,F5. 3, 4X, ’DEGREES’ ,/,lH  , 

+  17X, ’EXIT  ANGLE  =  ’ ,F8. 3 ,4X ,’ DEGREES ’,/, 1H  , 17X, 

+  ’ TURNNING  ANGLE  =  ’ ,F8. 3 ,4X ,’ DEGREES ’,/, 1H  ,17X, 
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+  'INLET  VEL  =  '  ,  F10.  6  , 3X , ' EXIT  VEL  =  ' ,F10. 6,//) 

WRITE 16, 36 5) 

365  FORMAT! 3X, ’ I ’ , 8X , ’ XO’ , 10X , ' YO ' , 10X, ' X ' , 11X, 'Y' , 10X, ' VT' , 

+  10X, ' VN1 ,11X, 'V' ,10X, 'CP' . 9X , ’ CPI ’ /// ) 

END  IF 

DO  360  1=1, N 

VTCOMi, I )=  VTAN( 1 , 1 )*COSAL+VTAN( 2 , 1 )*SINAL-!-W*VTANf  3,1) 

VNCOMC I  )=  VNORC 1 , 1  )*C0SAL+VN0R(.  2 , 1  )*SINAL+W*VNOR(  3,1) 

VC0M(  I)  =  SQRT(VTCOM(IV-">2  +  VNCOM(I)**2) 

IF  (.VTCOM(I)  .  LT.  0.0)  VCOM(I)  =  -VCOM(I) 

CP(I)  =  1.0  -  VCOM(I)  **  2 
CPI(I)=  1.0  -  ( VCOM( I )*UOUI )**2 

SIGCOM(I)  =  SIG(i,I)*COSAL+SIG(2,I)*SINAL+W*SIG(3,I) 

XPfl)  =0.5  *  (X0CI)+X0(I+1)) 

YP(I)  =0.5  *  ( YOC I )+YO( 1+1) ) 

380  CONTINUE 

IF  (ICYCLE  .  EQ.  1  .OR.  ICYCLE  .  GE.  ICYTL-1  -OR.  IP  .  GE.  0)  THEN 
WRITE  (1,374)  (  XO(I),YO(I),XP(I),YP(I),CP(I),CPI(I)  ,1=1, N) 

WRITE  (6,375)  (  I,  XO( I ) ,  YO( I ) ,  X?(I),  YP( I ) ,  VTCOM(I), 

+  VNCOM( I) ,VCOM( I) ,CP( I) ,  CPI(I)  ,I=1,N) 

3  74  FORMAT! 6F 10.  4) 

375  FORMAT! IX,  13,  9F12. 4) 

C  WRITE  (2)  I ,XO( I ) , YO( I ) ,XSTGR( I ) , YSTGR( I ) ,DLS( I),X(I),Y(I) , VC0M( 
WRITE! 6 , 385 )  N+i,XO(N+l) ,YO(N+l) 

335  FORMAT( IX, 13, 2F 12.  4) 

WRITE! 6 ,400)  CHORD,  CL1  ,CLI 

400  FORMAT! ///3X, ’CHORD  =  ’ ,F10.  5 ,4X, ’ CL(AVG)  =  ’ .F10.5.4X, 

+  ’CL( INLET)  =  1 , F10. 5) 

END  IF 

420  FORMAT! / 3X , 1HI , 6X , 2HSO , 10X , 2HSC , 9X , 3H  VN , 9X , 3HVNC , 9X , 3HDLS , 8X , 

+  4HDLSC) 

430  FORMAT! 15 , 6E12.  4) 

RETURN 

END 


DATA  SET  KCBCAMEAN  AT  LEVEL  001  AS  OF  08/24/84 
DATA  SET  KC3CAMEAN  AT  LEVEL  003  AS  OF  04/05/84 
SUBROUTINE  AMEAN( NS , ND , X , Y , IT) 

SMOOTH  DATA  USING  3-PTS  WEIGHTING  FORMULA 
NS  :  STARTING  PINT  OF  THE  DATA  TO  BE  SMOOTHED 

ND  :  END  PINT  OF  THE  DATA  TO  BE  SMOOTHED 

X,  Y  :  INDEPENDENT  +  DEPENDENT  VARAIBLES  OF  THE  DATA 
TO  BE  SMOOTHED 

IT  :  CYCLES  OF  DATA  SMOOTHING 


DIMENSION  X!101),Y!101) 


NM  =  ND  -NS 

I F ( NM  . LT.  2  .  OR.  IT  . LT. 1)  RETURN 

NDM 1  =  ND  -  1 
NSP1  =  NS  +  1 
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DO  20  K=i , IT 

DL1  =  Xi.NSPl)  -  X(NS) 

VI  -  Y(XS) 

DO  10  I=NSPi ,NDM1 
DL2  =  X( 1  +  n  -X(I) 

Y2  =  Y( I) 

YM  =  ( DL2  '  Y1  +  DL1  *  Y(I+1))/(DL1  +  DL2) 

YC I )  =0.5  *  ( Y2  +  YM) 

DL1  =  DL2 

Y1  =  Y2 

10  CONTINUE 

20  CONTINUE 

C 

RETURN 

END 

C  DATA  SET  KCBCBLGRID  AT  LEVEL  001  AS  OF  08/24/84 

C  DATA  SET  KCBCBLGRID  AT  LEVEL  001  AS  OF  08/24/84 

C  DATA  SET  KCBCBLGRID  AT  LEVEL  004  AS  OF  04/05/84 

SUBROUTINE  3LGRID(N,X,T,D1) 

C 

C  GENERATE  B.  L.  X-WISE  GRID  USING  MODIFIED  COSINE  DISTRIBUTION 
C 

DIMENSION  X( 101) ,T( 101) ,D1( 101) 

DATA  CRAD/57. 2957795/,  BPI/3. 14159265/ 

C 

C  . 

c 

NN  =  2  *  N  -  1 

EN  =  FLOAT((NN-l)/2) 

THO  =  10.  /CRAD 

CTOl  =  1.  +  COS(THO) 

DTH  =  (BPI  -  THO)  /  EN 

FI  =  FLOAT(N  -  2) 

DO  10  I=N ,NN 
FI  =  1.  0  +  FI 

II  =  I  -  N  +  1 

XII  =  THO  +  FI  *  DTH 

X( II)  =  (1.  0  +  COS(XII))/CT01 
10  CONTINUE 

XI  =  X( 1) 

XN  =  X(N) 

CH  =  XN  -XI 

FN1  =  FLOAT(N-l) 

N10  =  N/ 10 

DO  20  1=1, N 

T(  I)  =  FL0AT(I-1)/FN1 

X( I)  =  (X( I) -X1)/CH 

20  CONTINUE 

C  CALL  SMFIT(N10 ,N,T,X,D1,N10) 

C  IF(X(2).  LT.  0.  35  *  X( 3) )  X(2)  =  0.  35  *  X(3) 

C  CALL  SMFIT( 1,N,T,X,D1,2) 

CALL  AMEAN( 1,N,T,X,N10) 

C 

RETURN 

END 

C  DATA  SET  KCBCBL2D  AT  LEVEL  001  AS  OF  08/24/84 
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DATA  SET  KCBCBL2D  AT  LEVEL  001  AS  OF  08/24/84 

DATA  SET  KCBCBL2D  AT  LEVEL  012  AS  OF  04/06/84 

SUBROUTINE  EL2D  (  ITR , ISWPT , SURFID) 

PROGRAM  CALCULATES  VISCOUS/' INVISCID  INTERACTION'  USING  HILBERT 
INTEGRAL. 


COMMON  /BLCO/  NX , NXT , NP , NPT , NTR , IT, INVRS ,NS , IP 

COMMON  /BLC1/  F( 101 ,2) ,U( 101,2) ,V( 101,2) ,W( 101,2) ,B( 101,2) 

COMMON  / ELC2/  DELF( 101) ,DELU( 101) ,DELV( 101) ,DELW( 101) 

COMMON  /BLC7/  C( 100, 100) ,D( 100) ,DB( 100) ,DBP( 100) ,UE0( 100) ,GI 
COMMON/ EDDY 1 /  RL , RX , SQRX , RXNTR , GMTR , GMTRS ( 100), ALFAS (100), 
l-  FFS(  100)  ,RTS(  100) ,  IEDY.NXSPT 

COMMON  /SMRY/  W( 100) , ITP( 100) , ISL( 100) ,DLS( 100) ,CF( 100) , 

I-  THT(  100)  ,NPSTR(  100) 

COMMON  /GTY  /  X( 101) ,UE( 100) ,P1( 100) ,P2( 100) ,CEL,CELH 
COMMON  /BONV/  ITMAX,EPSL,EPST,CONV 

COMMON  /SAVE/  FS( 101) ,US( 101) ,VS( 101) ,WS( 101) ,BS( 101) 

COMMON  /BLIN/  TITLE( 20) ,XC( 100) ,YC( 100) , ISG( 100) ,DELS( 100) , 

1-  XCTR ,  XTR ,  I STRP ,  I  CYC  LE ,  ICYTL ,  XCTRS  (2)  ,TRFIND(  2) 

COMMON  /I SURF/  ISF 
COMMON/PLOT /v”:v 2) ,NXVP(20,2) ,ICC 
DIMENSION  SURFID(4) 


GENERATE  B.  L.  GRIDS  +  SET  INITIAL  CONDITIONS 

DO  5  I  =  1  ,  NXT 
ALFAS ( I)  =  0.0 
FFS(I)  =  1.  0 
RTS(I)  =1.0 
CONTINUE 

CALL  INPUT( ITR, ISWPT, SURFID) 

CALCULATE  HILBERT  COEFFS. ,  C(I,J) 

CALL  CALCIJ(NXT,0) 

LOOP  OF  CALCULATIONS 


20 


30 


NSS  =  NS 
NXSPT  =  NXT  +  1 

IF  (  ICYCLE  .EQ.  1  )  NS  =  NXT  +  1 
NX  =  NX  +  1 

CEL  =0.5  *  (X(NX)  +  X(NX-l))  /(X(NX)  -X(NX-l)) 

CELH  =  0.5  *  CEL 
IT  =0 

RX  =  UE ( NX ) *X( NX ) *RL 
SQRX  =  SQRT(RX) 

IT  =  IT  +  1 

IF( IT  .LE.  ITMAX)  GO  TO  40 

NXM1  =  NX-1 

CALL  HEADER(  TITLE .SURFID , I STRP  ) 

WRITE(6,  170  )  (M,X(M) ,CF(M) ,DLS(M) ,UE(M) ,P2  (M),THT(M), 

+  D(M),ALFAS(M),  ITP(M) ,NPSTR(M) ,M=1 ,NXM1) 

WRITEC  6 ,  160  )  NX 


INT06650 
INT06660 
INT06670 
INT06680 
INT06690 
INT06700 
INT067 10 
INT06720 
1NT06730 
INT06740 
INT06750 
INT06760 
INT06770 
INT06780 
1NT06790 
INT06800 
INT06810 
INT06820 
INT06830 
INT06S40 
INT06850 
INT06860 
INT06870 
INT06880 
INT06890 
INT06900 
INT06910 
INT06920 
INT06930 
INI 06940 
INT06950 
INT06960 
INT06970 
INT06980 
INT06990 
INT07000 
INT07010 
INT07020 
INT07030 
INT07040 
INT07050 
INT07060 
INT07070 
INT07080 
INT07090 
INT07100 
INT07110 
INT07120 
INT07130 
INT07140 
INT07150 
INT07160 
INT07170 
INTO 7 180 
INT07190 
INTO 7 200 


76 


STOP 

40  CONTINUE 

IFC  NX  .  GT.  NTP.)  CALL  EDDY 

CALL  COEFTR 
CALL  SOLY3 

IF( V( 1,2). Gl. 0. 0)  GOTO  60 
C 

C  EXTRAPOLATE  CALCULATED  D  FOR  TURBULENT  SEPARATION  OR  LAMINAR 
C  SEPARATION  FOR  LAMINAR  FLOW  CALCULATION  ONLY 

C 

CALL  EXTRAP ( NX , NXT , X , D ) 

NXM1  =  NX  -  1 

CALL  HEADER(  TITLE , SURFID , ISTRP  ) 

WRITE ( 6 ,  170  )  (M,X(M) ,CF(M) ,DLS(M) ,UE(M) ,P2(M) ,THT(M) , 

+  D(M) , ALFAS(M) , ITP(M) ,NPSTR(M) ,M=1,NXM1) 

WRITE( 6 , 180) 

KRITE(6, 190)  (M,X(M) ,D(M) ,M=NX,NXT) 

GOTO  130 

60  I F ( NX  .GT.  NTR)  GO  TO  70 

IF(ABS(DELY( 1) )  . GT.  EPSL)  GO  TO  30 
GO  TO  80 
70  CONTINUE 

IF(ABS(DELV( 1)/V(1,2))  . GT.  EPST)  GO  TO  30 
80  CONTINUE 

C 

C  CHECK  FOR  GROWTH 

C 

IF(NP  .  GE.  NP'7)  GO  TO  90 

IF(ABSCV(NP,2; )  .  LT.  0.0005  .AND.  ABS( 1. 0-U(NP-2,2) ) 

+  . LT. 0. 0035)  GOTO  90 

CALL  FILLUP(l) 

IT  =  1 
GO  W  30 
90  CONTINUE 

C 

CALL  FILLUP( 2) 

CALL  OUTPUT(l) 

I F ( ITR. EQ. 0  . OR.  NX. GE.  NTR)  GOTO  100 
IFCNX.LT.  3  .OR.  ITR.NE.  3)  GO  TO  100 
CALL  TRNS( I CODE) 

I F ( ICODE. EQ. 1)  GOTO  20 
100  IFCNX  . NE.  NSS)  GOTO  120 
C 

C  STORE  PROFILES  AT  THE  STATION  NS  FOR  INVERSE  B.  L. 

C  CALCULATION 
C 

DO  110  J  =  1  ,  NPT 
FS(J)  =  F( J, 2) 

US(J)  =  U(J,2) 

VS( J)  =  V( J , 2) 

WS(J)  =  W( J, 2) 

BS(J)  =  B(J,2) 

110  CONTINUE 

120  IF  (  NX  .LT.  NSS  )  GOTO  10 

C  IF  (ICYCLE  .NE.  1)  GO  TO  130 

IF  (  NX  .GE.  NS  )  GOTO  130 


IXT07210 
INT07229 
INT07230 
INTO 7240 
INTO  7  250 
INT07260 
INT07270 
INT07280 
INTO 7 290 
INT07300 
INT07310 
INT07320 
INT07330 
INT07340 
INTO? 350 
INT07360 
INT07370 
INT07380 
INT07390 
INTO  7400 
INT07410 
INTO 7420 
INT07430 
INTO  7  440 
INTO 7450 
INTO 7460 
INT07470 
INTO 7480 
INTO 7490 
INTO  7  5  00 
INT075 10 
INT07520 
INT07530 
INT07540 
INTO 7 550 
INTO 7 560 
INT07570 
INT07580 
INT07590 
INTO 7 600 
INT07610 
INTO 7620 
INT07630 
INT07640 
INT0765O 
INTO 7 660 
INT07670 
INT07680 
INT07690 
INT07700 
INT07710 
INTO 7 720 
INT07730 
INT07740 
INTO 7 750 
INT07760 


77 


IF  (NX  .  LT.  NXT)  GO  TO  10  INT07770 

CALL  HEADER(  TITLE , SURFID , ISTRP  )  IST07780 

WRITE(6,  170  )  (M,X(M) ,CF(M) ,DLS(M) ,UE(M) ,P2  (M) ,THT(M) ,  INT07790 

+  DU!),  ALFAS(M) ,  ITP(M) ,NPSTR(M) ,M=1 ,NXT)  IST07800 

130  DO  140  I  =  1  ,  NXT  INT07810 

140  DB( I )  =  D( I )  INT07820 

NS  =  NSS  INT07830 

NX  =  NS  INT07840 

NP  =  NPSTR(NX)  INT07850 

DO  150  J  =  1  ,  NPT  INT07860 

F( J,2)  =  FS( J)  INTO 78 70 

U ( J , 2 )  =  US(J)  INT07880 

V( J , 2)  =  VS(J)  INT07890 

V ( J , 2 )  =  WS(J)  INTO 7900 

B(J,2)  =  BS( J)  INT07910 

150  CONTINUE  INT07920 

155  INVRS  =  NS  +  1  INTO 7 9 30 

C  INT07940 

C  CALCULATION  SHIFTS  TO  USING  PHYSICAL  COORDINATES  INT07950 

CALL  MATN2(ITR,ISVPT, SURFID)  INT07960 

C  INT07970 

C  PASS  DELTA-STAR  BACK  TO  MAIN  PROG.  INT07980 

C  INT07990 

DO  158  I  =  1 , NXT  INT08000 

DELS(I)  =  DLS(I)  INT08010 

158  CONTINUE  INT06020 

RETURN  INT08030 

C  INT08040 

C  .  INT08050 

C  INT08060 

160  FORMAT( 1H0 , '  **  ITERATIONS  EXCEEDED  ITMAX  AT  NX  =  ',15/  INT08070 

+  1H  **  CALCULATIONS  STOP.  **  ')  INT08080 

170  F0RMATC1H0,'  **  SUMMARY  OF  STANDARD  B.  L.  SOLUTIONS.  **'/  INT08090 

+  1H0 , 4X , 2HNX , 7X , 1HX , 12X , 2HCF , 1 IX , 3HDLS , 12X , 2HUE ,  INT08100 

+  12X, 2HP2 , 1 IX , 3HTHT, 13X, 1HD, 10X.4HALFA, 6X, 2HIT, 2X, 2HNP/IMT08110 

+  (1H  ,3X,I3,F10.  5,2X,7E14. 5,214))  INT08120 

180  FORMAT( 1H0 , 34H  FLOW  SEPARATES.  D  IS  EXTRAPOLATED/  INT08130 

+  1H0 , 3X, 3H  NX , 7X , 1HX , 1 3X , 1HD/ )  INT08140 

190  FORMAT( 1H  ,3X  13 ,F10.  5 , 2X,E14. 5)  INT08150 

END  INT08160 

C  DATA  SET  KCBCCALCIJ  AT  LEVEL  001  AS  OF  08/24/84  INT08170 

C  DATA  SET  KCBCCALCIJ  AT  LEVEL  001  AS  OF  08/24/84  INT08180 

C  DATA  CE7  KCBCCALCIJ  AT  LEVEL  005  AS  OF  04/05/84  INT08190 

SUBROUTINE  CALCIJ  (  IL,  LO)  INT08200 

C  INT08210 

C  CALCULATE  HILBERT  INTEGRAL  COEFFS  INT08220 

C  INT08230 

COMMON  /BLC7/  C( 100 , 100) ,D( 100)  ,DB( 100)  ,DBP( 100) ,UEO( 100) ,GI  INT08240 

COMMON/EDDY 1/  RL,RX,SQRX,RXNTR,GMTR,GMTRS( 100)  INT08250 

+  , ALFAS( 100) ,FFS( 100) ,RTS( 100) , IEDY ,NXSPT  1NT08260 

COMMON  /GTY  /  X( 101) ,UE( 100) ,P1( 100) ,P2( 100) ,CEL,CELH  INT08270 

DIMENSION  E( 2)  INT08280 

C  INT08290 

C  .  INT08300 

C  INT08310 

PI  =  3. 14159265  INT0S320 
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c 

c 

c 


c 

c 

c 

30 


c 

c 

c 

40 

C 

50 


60 


65 

C 

c 


c 

c 

c 

c 

c 

c 

c 


PI  =  PI*SQRT(RL) 

IL1  =  IL  -  1 
DO  65  I  =  2,  IL1 
E  (1)  =  0. 

L  =  LO  +  I 


DO  60  J  = 
J1 
K 

DX1 
DX2 
DX3 
IF  ( 

IF  ( 


o 

<-  > 


J 

J 


IL 
=  J  -  1 
=  J  +  LO 
=  X(L)  - 
=  X(K)  - 
=  X(L)  - 
EQ. 

EQ. 


X(K) 
X(K-l) 
X(K-l) 

I  )  GO  TO  30 
(1+1)  )  GO  TO 


40 


J  . NE.  I  OR  1+1 


E  (2) 
GO  TO 


=  (  1.0/DX2  )  *  ALOG(  ABS(  DX3  /  DX1  )  ) 


50 


J  .EQ.  I 


R1 

E  (2) 
GO  TO 


X(K+1) -X(K)  )  / 
R1  *  ALOG(  ABS( 


(X(K+1)-X(K-1)  ) 

DX3  /  (  X(L)-X(K+1)) 


50 


J  .EQ.  1+1 

R1  =  ( 

E  (2)  =  < 

CONTINUE 
C  (J1,I)=  l 
E(  1 ) 
CONTINUE 
E  (2)  = 

J1 

K  = 

C  (J1,I)= 
CONTINUE 


X(K 
R1  • 


•l)-X(K-2)  ) 
■  ALOG(  ABS( 


/  (  X(K)-X(K-2) 
(X(L) -X(K-2) )  / 


) 

DX1 


=  £, 


E(  1) 

(2) 


-  E(2)  )  /  PI 


0. 

IL 

K  +  1 
E(  1)  / 


PI 


RETURN 

END 

DATA 

DATA 

DATA 

SUBROUTINE 


SET 

SET 

SET 


KCBCCOEF 

KCBCCOEF 

KCBCCOEF 


AT 

AT 

AT 


LEVEL 

LEVEL 

LEVEL 


001 

001 

007 


AS 

AS 

AS 


OF 

OF 

OF 


COEF( GAMMA 1 , GAMMA2) 


CALCULATE  COEFFS  OF  B. 
SEMI  -TRANSF  VAR.' TBLES( 


L.  FINITE-DIFFERENCE  EQS. 
AFTER  SWITCHING). 


COMMON  /BLCO/  NX , NXT , NP , NPT , NTR , IT , INVRS , NS , IP 
COMMON  /BLC1/  " 1101 , 2) ,U( 101 , 2) , V( 101 , 2) ,V( 101 , 2) , B( 101 , 2) 

COMMON  /ELU6/  Sl( 101) ,S2( 101) ,S3( 101) ,S4( 101) ,S5( 101) ,S6( 101) , 
h  S7(101),S8(101),R1(101),R2(101),R3(101),R4( 101) INT0S870 

COMMON  /BLC7/  C( 100 , 100) ,D( 100) ,DB( 100) ,DBP( 100) ,UEO( 100) ,GI  INT08880 


INT08330 

INT08340 

INT08350 

INT08360 

INT08370 

INT08380 

INT08390 

INT08400 

INT08410 

INT08420 

INT08430 

INT08440 

INT08450 

INT08460 

INT08470 

INT08480 

INT08490 

INT08500 

INT08510 

INT08520 

INT08530 

INT08540 

)  )  +  2. 0  )  /  DX2INT08550 
INT08560 
INT06570 
INT08580 
INT08590 
INT08600 

)  )  -  2.  0  )  /  DX2  INT08610 
INT08620 
INT08630 
INT08640 
INT08650 
INT08660 
INT08670 
INT08680 
INT08690 
INT08700 
INT087 10 
INT08720 
INT08730 
INT08740 
INT08750 
INT08760 
INT08770 
INT08780 
INT08790 
INT08800 
INT08810 
INT08820 
INT08830 
INT08840 
INT08850 
INT08860 


08/24/84 

08/24/84 

04/05/84 


IN 


79 
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COMMON  /GRD  /  ETA( 101) ,DETA( 101) ,A( 101) 

COMMON  /GTY  /  X( 101) ,UE( 100) ,P1( 100) ,P2( 100) ,CEL,CELH 


P1H  =0.5  *  PI (NX) 

DO  100  J=  2 ,NP 
FLARE  =1.0 

FB  =  0.  5*(F( J,2)  +  F( J- 1 , 2) ) 

UB  =  0. 5*(U( J,2)  +  U( J- 1,2)) 

FVB  =  0. 5*(F(J,2)*V(J,2)+F(J-1,2)*V(J-1,2)) 

IF(UB  .LT.  0.0)  FLARE  =  0.0 
VB  =  0.  5*(V(J,2)  +  V( J-1,2)) 

USB  =  0.  5*(U(J,2)**2  +  U( J-l,2)**2) 

WSB  =  0. 5*(W(J,2)**2  +  W( J-1,2 )**2) 

DERBV  =(B(J,2)*V(J,2)  -  B( J- 1 , 2)*V( J-l , 2) )/DETA( J-l) 

FB4  =  0. 5*(F( J , 1)  +  F( J- 1 , 1)  ) 

VB4  =  0.  5*(V(J,1)  +  V(J-1,1)) 

USB4  =  0.  5*(U(J,1)**2  +  U( J-l , 1)**2) 

V.SB4  =  0.  5*(W(J,1)**2  +  W(J-1,1)**2) 

FVB 4  =  0.5*(F(J,lj*V(J,l)+F(J-l,l)*V(J-l,l)) 

DERBV4  =  (B( J, 1)*V( J, 1)  -  B( J-l , 1)*V( J-l , 1) )/DETA( J-l) 

S1(J)  =  CELH*(FB  -  FB4)  +  P1H*F(J,2)  +  B( J,2)/DETA( J-l) 

S2( J)  =  CELH*(FB  -  FB4)  +  P1H*F(J-1,2)  -  B( J-1,2)/DETA( J-l) 

S3(J)  =  CELH*( VB  +  VB4)  +  P1H*V(J,2) 

S4( J)  =  CELH*( VB  +  VB4)  +  P1H*V(J-1,2) 

S5( J)  =  -CEL*FLARE*U( J , 2 ) 

S6( J)  =  -CEL*FLARE*U(J-1,2) 

S7( J)  =  CEL" W ( J , 2 ) 

S8( J)  =  CEL*W( J-l ,2) 

CRB  =  -DERBV4  +  CEL*WSB4  -  CEL*FLARE*USB4  -  P1(NX)*FVB4 
R2(J)  =  CRB  -  (DERBV  -  CEL*FLARE*USB  +  CEL*(VB+VB4)*(FB-FB4)  + 

+  CEL*WSB  +  P1(NX)*FVB) 

R1(J)  =  F( J-1,2)  -  F(J,2)  +  DETA(J-1)*UB 

R3( J-l)=  U( J-1,2)  -  U(J,2)  +  DETA(J-1)*VB 
R4( J-l)=  W( J-1,2)  -  W ( J , 2 ) 

CONTINUE 

BOUNDARY  CONDITIONS 

Rl( 1)  =  0. 0 

R2( 1)  =  0. 0 

R4(NP)  =  0. 0 

IF(NX  . GE.  INVRS)  GO  TO  120 
GAMMA  1  =  0.0 
GAMMA2  =  1.0 
R3(NP)  =  0. 0 
RETURN 
120  CONTINUE 

CII  =  C(NX,NX)  *  SQRT(X(NX)) 

GAMMA  1  =  1.0 

GAMMA 2  =  (1.0  -  CII*ETA(NP) )/CII 

R3( NP)  =  (GI  +  C I I - ( ETA( NP ) *W ( NP , 2 )  -  F(NP,2))  -W(NP,2) )/CII 
C 

RETURN 


INT08890 
INT08900 
INT08910 
INT08920 
INT08930 
INT08940 
INTO 89 50 
INT08960 
INT08970 
INT08980 
INT08990 
INT09000 
INT09010 
INT09020 
INT09030 
INT09040 
INT09050 
INT09060 
INT09070 
INT09080 
INT09090 
INT09100 
INT09110 
INT09120 
INT09130 
INT09140 
INT09150 
INT09160 
INT09170 
INT09180 
INT09190 
INT09200 
INT09210 
INT09220 
INT09230 
INT09240 
INT09250 
INT09260 
INT09270 
INT09280 
INT09290 
INTO 9 3 00 
INT09310 
INT09320 
INT09330 
INT09340 
INT09350 
INT09360 
INT09370 
INT09380 
INT09390 
INT09400 
INT09410 
INT09420 
INT09430 
INT09440 


SO 
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DATA  SET  KCBCCOEFTR  AT  LEVEL  001  AS  OF  08/24/84 

DATA  SET  KCBCCOEFTR  AT  LEVEL  001  AS  OF  08/24/84 

DATA  SET  KCBCCOEFTR  AT  LEVEL  004  AS  OF  02/21/84 

SUBROUTINE  COEFTR 

CALCULATE  COEFFS.  OF  B.  L.  FINITE-DIFFERENCE  EQS. 

IN  TRANSFORMED  VARIABLES(  BEFORE  SWITCHING). 

COMMON  /BLCO/  NX , NXT , NP , NPT , NTR , IT , INVRS , NS , IP 
COMMON  /BLC1/  F( 101,2) ,U( 101,2) ,V( 101,2) ,W( 101,2) ,B( 101,2) 
COMMON  /BLC6/  Sl( 101) ,S2( 101) ,S3( 101) ,S4( 101) ,S5( 101) ,S6( 101) , 
+  S  7 ( 101) ,SB(101) ,R1( 101) ,R2( 101) ,R3( 101), R4( 

COMMON  /GRD  /  ETA( 101) ,DETA( 101) ,A( 101) 

COMMON  /GTY  /  X( 101) ,UE( 100) ,P1( 100) ,P2( ICO) ,CEL,CELH 


DO  100  J=  2 ,NP 

FB  =  0.5*(F(J,2)  +  F( J- 1 , 2) ) 

UB  =  0. 5*(U(J,2)  +  U( J-1,2)) 

VB  =  0. 5*(V(J,2)  +  V( J-1,2)) 

USB  =  0. 5*(U(J,2)**2  +  U( J-l ,2)**2) 

DERBV  =  ( B ( J , 2)*V( J , 2)  -  B( J-l , 2)*V( J-l , 2) ) /DETA( J-l) 
FVB  =  0. 5*(F(J,2)*V(J,2)  +  F(J-1,2)*V( J-1,2)) 

FV’34  =  0.  5*(F(J,1)*V(J,1)  +  F(  J-l ,  1)*V( J-l ,  1)) 

FB4  =  0. 5*(F( J, 1)  +  F( J-l , 1) ) 

VB4  =  0. 5*( V( J, 1)  +  V(J-1,1)) 

USB4  =  0. 5*(U(J,1)**2  +  U( J-l , 1)**2) 

DERBV4  =  ( B ( J , 1 ) *V( J , 1 )  -  B( J-l , 1)*V( J-l , 1) )/DETA( J-l) 


S1(J)  =  CELH*(FB-FB4)  +  0. 5*P1 
S2( J)  =  CELH*(FB-FB4)  +  0. 5*P1 
i-  DETA(J-l) 

S3( J)  =  CELH*(VB  +  VB4)  +  0.5* 
S4( J)  =  CELH*(VB  +  VB4)  +  0.5* 
S5( J)  =  -(CEL+P2(NX))*U(J,2) 
S6( J)  =  -( CEL+P2( NX) )*U( J-1,2) 


0.  5*P1(NX)*F(J,2)  +  B( J,2)/DETA( J-l) 
0.  5*P1(NX)*F(  J-1,2)  -B( J-l ,2)/ 

+  0. 5*P1(NX)*V(J,2) 

+  0. 5*P1(NX)*V( J- 1,2) 


CLB  =  DERBV4  +  P1(NX-1)*FVB4  -  P2(NX-1)*USB4  +  P2(NX-1) 

CRB  =  -CLB  -  CEL*USB4  -  P2(NX) 

R2(J)  =  CRB  -  (DERBV  +  P1(NX)*FVB-  ( CEL+P2(NX) )*USB  +  CEL* 
1-  (VB  +  VB4)*(FB  -  FB4) ) 

R1(J)=  F( J-1,2)  -  F( J , 2)  +  DETA( J-1)*UB 

R3( J-l)=  U( J-1,2)  -  U( J , 2)  +  DETA( J-1)*VB 

CONTINUE 

Rl(  1 )  =  0.  0 

R2(  1)  =  0.0 

R3(NP)  =0.0 

RETURN 

END 

DATA  SET  KCBCCOMPBL  AT  LEVEL  001  AS  OF  08/24/84 
DATA  SET  KCBCCOMPBL  AT  LEVEL  001  AS  OF  08/24/84 


INT09450 
INT09460 
INT09470 
INT09480 
INT09490 
INT09500 
INT095 10 
INT09320 
INT09530 
INT09540 
INT09550 
INT09560 
101) INT09570 
INT09580 
INT09590 
INT09600 
INT09610 
INT09620 
INT09630 
INT09640 
INT09650 
INT09660 
INT09670 
INT09680 
INT09690 
INT09700 
INT09710 
INT09720 
INT09730 
INT09740 
INT09750 
INT09760 
INT09770 
INT09780 
INT09790 
INT09800 
INT09810 
INT09820 
INT09830 
INT09840 
INT09S50 
INT09860 
INT09870 
INT09880 
INT09890 
INT09900 
INT09910 
INT09920 
INT09930 
INT09940 
INT09950 
INT09960 
INT09970 
INT09980 
INT09990 


81 


c 

r> 


C 


c 


c 


c 

c  - 

c 

90 

110 

112 

120 

122 

130 

140 


150 

160 

170 

180 

182 

190 

200 

220 

210 

C 

C 

C 


DATA  SET  KCBCCOMPBL  AT  LEVEL  010  AS  OF  08/24/84 


SUBROUTINE  COMPBL( CASE ID , XP , YP , VT , S , DLSP , DLS , DBPP , NBL , ITRI , XCTRI 
+  RN,NT, ISWPT) 


COMMON  /BLCO/  NX ,NXT,NP,NPT,NTR, IT, INVRS ,NS , IP 

COMMON  / B LC 7 /  C( 100 , 100) ,D( 100) ,DB( 100) ,DBP( 100) ,UE0( 100) ,GI 

COMMON/ EDDY 1/  RL , RX , SQRX , RXNTR , GMTR , GMTRS ( 100) 

+  , ALFAS( 100) ,FFS( 100) ,RTS( 100) , IEDY.NXSPT 

COMMON  /GTY  /  X( 101) ,UE( 100) ,P1( 100) ,P2( 100) ,CEL,CELH 
COMMON  /BLIN/  TITLE( 20) ,XC( 100) ,YC( 100) , ISG( 100) ,DELS( 100) , 

+  XCTR,XTR, ISTRP , ICYCLE , ICYTL,XCTRS( 2) ,TRFIND( 2) 

COMMON  /BLOW/  VN(IOO) 

COMMON  / I SURF/  ISF 
COMMON/PLOT/NVP( 2) ,NXVP( 20,2) ,ICC 


DIMENSION  XP( 100) ,DLSP( 100) ,YP( 100) ,VT( 100) ,S( 100) , 
t-  DBPP(  100)  ,DLS(  100) ,  CASEID(  20) 

DIMENSION  XIN( 100,2) ,YIN( 100,2) ,SI( 100,2) ,VIN( 100,2) , 
t-  DIN(  100,2)  ,DELSTR(  100,2)  ,DD(  100,2)  ,DDD(  100,2) 

DIMENSION  XB( 101) ,D1(101) ,D2(101) ,D3( 101) ,IEND(2) 

DIMENSION  NBL( 2) , ITRI( 2) ,XCTRI(2) 

LOGICAL  TRFIND 

CHARACTER  *  4  SURF(4) ,STITLE(2) ,SURFID(4) 


DATA  SURF  /  ’  ’ , 'R  SU’ , 'RFAC' , 'E  '  / 

DATA  STITLE  /  ' UPPE ' , ' LOWE '  / 


FORMAT  (  1H1,5X, 'COMPUTING  BOUNDARY  LAYER  USING  HILBERT' , 

+  '  INTEGRAL',/) 

FORMAT  (  1H0,6X,'I',9X,'XP',13X,'YP',13X,'S',14X,'VT',13X, 

+  ' DBP'  /  ) 

FORMAT  (  1H0 , 6X, ' I ' ,4X, ' II ' , 3X, ' IK' , 7X, ' XIN' , 12X, ' YIN1 , 

+  13X, ' SI ' , 12X, ' VIN' , 12X, 'DIN'  /  ) 

FORMAT  (  1H  ,5X,I3,5E15. 6  ) 

FORMAT  (  1H  , 3X, 3( 2X, 13) , 5E15. 6  ) 

FORMAT  (  1H0,5X, 'STAGNATION  POINT  IS  FOUND  BETWEEN  POINT  NO.  ', 
+  13, '  AND  POINT  NO.  ’ ,13  /  ) 

FORMAT  (  1H0,5X, ' INTERPOLATED  STAGNATION  POINT  VALUES'  / 

+  1H0 ,5X, ' S  =  ’ ,E13. 6,2X,'XP  =  ' ,E13.  6,2X, 'YP  =  ', 

+  E13. 6,2X, 'DBP  =  ’ ,E13.  6,2X, ’VT  =  0.0’  /  ) 

FORMAT  (  1H0,5X, 'TOTAL  NUMBER  OF  UPPER  SURFACE  POINTS  =  ',15, 

+  2X , ' AND  AT  LOWER  SURFACE  =  ',15  /  ) 

FORMAT  (  1H0.5X, 'UPPER  SURFACE  DATA'  ) 

FORMAT  (  1H0,5X, 'LOWER  SURFACE  DATA'  ) 

FORMAT  (  1H0,5X, 'UPPER  SURFACE  CALCULATIONS'  /  ) 

FORMAT  (  1H0,5X, 'LOWER  SURFACE  CALCULATIONS'  /  ) 

FORMAT  (  1H0,5X, 'RESULTS  OF  POINT  REDISTRIBUTION'  /  ) 

FORMAT  (  1H0,5X, 'TABLE  OF  DELTA-STARS’  /  (  1H  ,5X,8E15.6  )  ) 
FORMAT  (1H0,5X, ’TABLE  OF  BLOWING-VEL. ’  /  ( 1H  ,5X,8E15.6)  ) 
FORMAT  (  1H0 , 5X, 1  NO  CHANGE  OF  SIGN  ON  VT.  CANNOT  FIND  STAG.  PT.  ' 


INT10000 
INT10010 
INT10020 
INT10030 
INT10040 
INT10050 
INT10060 
INT10070 
INT10080 
INT10090 
INT10100 
INT10110 
INT10120 
INT10130 
INT10140 
INT10150 
INT10160 
INT10170 
INT10180 
INT10190 
INT10200 
INT10210 
INT10220 
INT10230 
I NT 10 240 
INT10250 
INT10260 
INT10270 
INT10280 
INT10290 
INT10300 
INT10310 
INT10320 
INT10330 
INT10340 
INT10350 
INT10360 
INT10370 
INT10380 
INT10390 
INT10400 
INT10410 
INT10420 
INT10430 
INT10440 
INT10450 
INT10460 
INT10470 
INT10480 
INT10490 
INT10500 
INT10510 
) INT10520 
INT10530 
INT10540 
INT10550 


S2 


c 

INT10560 

c 

READ  ONE  STRIP  INPUT  DATA  FROM  UNIT  NO.  1.  THE  ORDER  IS 

INT10570 

c 

FROM  THE  LOWER  SURFACE  T.  E.  TO  THE  UPPER  SURFACE  T. E. 

INT10580 

c 

INT10590 

DO  230  I  =  1,20 

INT10600 

TITLE(I)  =  CASEID(I) 

INTI 06 10 

230 

CONTINUE 

INT10620 

RL  =  RN 

INT10630 

DO  300  I  =  1,NT 

INT10640 

DBP(I)  =  DBPPvT) 

INT10650 

300 

CONTINUE 

INT10660 

C 

INT10670 

C 

PRINT  THE  INPUT  DATA. 

INT10680 

C 

INT10690 

P2(n  =  1.0 

INT10700 

C 

WRITE  (  6,90  ) 

INT107 10 

C 

WRITE  (  6,110  ) 

INT10720 

c 

DO  500  I  =  1 ,NT 

INT10730 

c 

WRITE  (  6,120  )  I,XP(n,YP(I),S(I),VT(I),DBP(I) 

INT10740' 

C500 

CONTINUE 

INTI 0750 

C 

INT10760 

C 

FIND  STAGNATION  POINT 

INT10770 

C 

INT10780 

DO  600  I  =  2, NT 

INT10790 

VPROD  =  VT( I )  *  VT(I-l) 

INT10800 

IF  (  VPROD  . GT.  0.  )  GO  TO  600 

INT10810 

IS  =  I 

INT10820 

IS1  =  IS  -  1 

INT10830 

GO  TO  700 

INT10840 

600 

CONTINUE 

INT10850 

WRITE  (  6,210  ) 

INT10860 

STOP  1 

INT10870 

C 

INT10880 

C 

INT10890 

C 

INTERPOLATE  S  AT  VT  =  0. 

INT10900 

C 

INT10910 

C 

WRITE  (  6,130  )  ISl.IS 

INT10920 

700 

DS  =  S(IS)  -  S(IS1) 

INT10930 

DV  =  VT(IS)  -  VT(ISl) 

INT10940 

SS  =  S(IS)  -  VT(IS)  *  (  DS/DV  ) 

INT10950 

DBB  =  DBP(IS)  -  DBP(ISl) 

INT10960 

DX  =  XP(IS)  -  XP(ISl) 

INT10970 

DY  =  YP(IS)  -  YP(ISl) 

INT1098U 

DS1  =  SS  -  S(IS) 

INT10990 

DBS  =  DBP(IS)  +  DS1  *  (  DBB/DS  ) 

INTI 1000 

XPS  =  XP(IS)  +  DS1  *  (  DX  /DS  ) 

1NT11010 

YPS  =  YP(IS)  +  DS1  *  (  DY  /DS  ) 

INTI 1020 

c 

WRITE  (  6,140  )  SS, XPS, YPS, DBS 

INTI 1030 

c 

INT11040 

c 

IU  IS  THE  TOTAL  UPPER  SURFACE  POINTS.  +  STAG.  PT. 

INTI  1050 

c 

IL  IS  THE  TOTAL  LOWER  SURFACE  POINTS  +  STAG.  PT. 

INTI 1060 

c 

INT11070 

IU  =  NT  -  IS  +  2 

INTI 1080 

IL  =  IS 

INT11090 

c 

WRITE  (  6,150  )  IU, IL 

INTI 1100 

c 

INTI 11 10 

83 


C  GROUP  TEE  DATA  FOR  EACH  SURFACE.  FIRST,  UPPER. 

C 

DO  1200  L  =  1,2 
GO  TO  (  800,900  ),  L 
C 

C  L  =  1  IS  UPPER  SURFACE 

C  L  =  2  IS  LOWER  SURFACE 

C 

800  Ml  =  IS 

M2  =  NT 

IEND(L)  =  IU 
GO  TO  1000 
900  Ml  =  1 

M2  =  IL-1 

IEND(L)  =  IL 
C 

1000  1=1 

XIN'C  I ,  L)  =  XPS 

YIN(I.L)  =  YPS 
SI  ( I , L)  =  0. 

DIN( I , L)  -  DBS 
VINCI.L)  =  0. 

IF  (  IP  .GE.  1  )  THEN 

IF  (  L  .EQ.  1  )  WRITE  (  6,160  ) 

IF  (  L  .EQ.  2  )  WRITE  (  6,170  ) 

WRITE  (  6,112  ) 

WRITE  (  6,122  )  I , I , I ,XIN( 1 ,L) ,YIN( 1 ,L) , SI( 1 ,L) , VIN( 1 ,L) , 

+  DIN( 1 , L) 

END  IF 

DO  1100  II  =  Ml, M2 
I  =1  +  1 

IK  =  II 

IF  (  L  .EQ.  2  )  IK  =  IL  -  II 
XIN( I ,L)  =  XP(IK) 

YIN( I , L)  =  YP(IK) 

SI  ( I ,L)  =  S( IK) -SS 

IF  (  L  .EQ.  2  )  SI( I ,L)  =  SS-S(IK) 

VIN( I ,L)  =  ABS(VT( IK) ) 

DIN( I , L)  =  DBP(IK) 

IF C IP  .GE.  1) WRITE  (  6,122  )  I , II , IK,XIN( I ,L) , YIN( I ,L) , SI( I ,L) , 
+  VIN(I,L),DIN(I,L) 

1100  CONTINUE 
C 

1200  CONTINUE 
C 

C  RE-DISTRIBUTE  POINTS  ON  EACH  SURFACE  TO  A  DENSER  NUMBER. 

C 

C  WRITE  (  6,90  ) 

DO  2000  ISF  =  1,2 
NN  =  IEND(ISF) 

ITR  =  ITRI(ISF) 

NXT  =  NBL(ISF) 

XCTR  =  XCTRI(ISF) 

SURF  (1)  =  STITLE(ISF) 

ICC  =  1 
DO  1220  J  =  1,4 


INTI 1120 
INT11130 
INTI 1140 
INTI 1 150 
INTI 1160 
INTI 11 70 
INT11180 
INTI  1 i90 
INT11200 
INT11210 
INT11220 
INTI 1230 
INTI  1240 
INT11250 
INT11260 
INT11270 
INTI  1280 
INTI  1290 
INTI 1300 
INT11310 
INTI  1320 
INT11330 
INTI 1340 
INTI  1350 
INTI  1360 
INT11370 
INT11380 
INT11390 
INTI  1400 
INT11410 
INT11420 
INTI  1430 
INTI  1440 
INT11450 
INT11460 
INT11470 
INT11480 
INTH490 
INTI  1500 
INT11510 
INT11520 
INT11530 
INTI  1540 
INTI 1550 
INT11560 
INT11570 
INT11580 
INT11590 
INTI  1600 
INT11610 
INTI 1620 
INT11630 
INTI 1640 
INTI 1650 
INT11660 
INT11670 


S4 


1220 

C 

C 

C 


n 

U 

1300 

C 

C 

c 


1320 

C 

C 

C 


1350 

C 

C 


c 

2000 

C 

C 

C 


SURFID(J)  =  SURF(J) 

CONTINUE 

IF  (  ISF  .EQ.  1  )  WRITE  (  6,180  ) 
IF  (  ISF  .EQ.  2  )  WRITE  (  6,182  ) 
SCALE  =  SI(NN, ISF) 

CALL  BLGRID  (  NXT,XB,D1,D2  ) 


DO  1300  I  =  1 , NXT 
X  (I)  =  XB( I)  *  SCALE 

INTERPOLATE  S,VT,X,Y,D  INTO  THE  NEW  DISTRIBUTION. 


CALL  SMFIT  (  1,NN,SI( 
CALL  SMFIT  (  1,NN,SI( 
CALL  DIFF3  (  NN,SI(1, 
CALL  INTRP3(  NN,SI(1, 
CALL  DIFF3  (  NN,SI(1, 
CALL  INTRP3(  NN,SI(1, 
CALL  DIFF3  (  NN,SI(1, 
CALL  INTRP3(  NN,SI(1, 
CALL  DIFF3  (  NN,SI(1, 
CALL  INTRP3(  NN,SI(1, 
IF  (  IP  .GE.  1  )  THEN 
WRITE  (  6,190  ) 

WRITE  (  6,110  ) 

DO  1320  I  =  1 .NXT 
WRITE  (  6,120  )  I,XC( 
CONTINUE 
END  IF 


1, ISF), VIN(1, ISF), Dl, 2  ) 

1 ,  ISF) ,DIN( 1 , ISF) ,D1 , 2  ) 

ISF), VIN(1, ISF), D1,D2,D3,0  ) 

ISF), VINCI, ISF) ,D1 ,D2 ,D3 ,NXT,X,UE  ) 
ISF) ,DIN( 1 , ISF) ,D1 ,D2 ,D3 , 0  ) 

ISF) ,DIN( 1 , ISF) ,D1,D2 ,D3,NXT,X,DB  ) 
ISF) ,XINC 1, ISF) ,D1 ,D2 ,D3 ,0  ) 

ISF) ,XIN( 1, ISF) ,D1 ,D2 ,D3 ,NXT,X,XC  ) 
ISF)  ,  YIN'C  1 ,  ISF)  ,D1  ,D2  ,D3,0  ) 

ISF) ,YINC 1,ISF) ,D1,D2,D3,NXT,X,YC  ) 


I),YCCI),XCI),UECI),DBCI) 


INPUT  TO  THE  B.  L.  PROGRAM  X,UE,DB,XC,YC  ARE  NOW  DEFINED. 

DO  1350  I  =  1 ,NXT 
DBPCI)  =  DB( I ) 

D( I )  =  DBCI) 

CONTINUE 


CALL  BL2D(  ITR, ISWFT, SURFID) 

CALL  DIFF3  (  NXT, X, DELS ,D1 ,D2 ,D3 , 0  ) 

CALL  INTRP3C  NXT ,X, DELS ,D1 ,D2 ,D3 ,NN, SI( 1 , ISF) ,DELSTR( 1 , ISF) 
CALL  DIFF3CNXT,X,D,D1,D2,D3,0) 

CALL  INTRP3CNXT,X,D,D1,D2,D3,NN,SIC1,ISF),DDC1,ISF)) 

CALL  DIFF3(NN, SIC  1 , ISF) ,DDC 1 , ISF) ,DDDC 1 , ISF) ,D2 ,D3 , 0) 
TRFINDC ISF)  =  .FALSE. 

I F  ( ITR  .EQ.  3  .AND.  NTR  .  LE.  NXT)  THEN 
XCTRS(ISF)  =  XCTR 
TRFINDC ISF)  =  .TRUE. 

END  IF 


CONTINUE 

PUT  THE  TWO  SURFACES  DELTA-STARS  BACK  TO  ONE  STRIP 
DELSTRC 1 ,2)  =  0.  5*( DELSTR( 2 , 1)+DELSTR( 2 , 2) ) 


INTI 1680 
INTI 1690 
INTI 1700 
INTI 17 10 
INTI 1720 
INTI  1730 
INT11740 
INTI 1750 
INT11760 
INTI 17 70 
INT11780 
INT11790 
INT11800 
INTI 18 10 
INTI 1820 
INT11830 
INT11840 
INTI 1850 
INTI 1860 
INT11870 
1NT11880 
INT11890 
INT11900 
INT11910 
INT11920 
INTI 1930 
INT11940 
INT11950 
JNT11960 
INT11970 
INTI 1980 
INT11990 
iNf i2000 
INT12010 
INT12020 
INT12030 
INT12040 
INT12050 
INT12060 
INT12070 
INT12080 
INT12090 
INT12100 
INT12110 
INT12120 
INT12130 
INT12140 
INT12150 
INT12160 
INT12170 
INT12180 
INT12190 
INT12200 
INT12210 
INT12220 
INT12230 


85 


non  non  ono 


DELSTR(l.l)  =  DELSTR( 1,2) 

DD( 1,2)  =  0.5  *  ( DD( 2,1)  +  DD(2,2)) 

DD( 11)  =  DD( 1,2) 

DDD( 1,2)  =  0.5  *  ( DDD( 2,1)  +  DDD(2,2))  /SQRT(RL) 

DDD( 1,1)  =  DDD( 1,2) 

IL  =  IEND( 2) 

I  =  IL 

DO  2100  II  =  2 , IL 
I  =  1-1 

DLS(I)  =  DELSTR( 11,2) 

VN( I )  =  DDD( 11,2) /SQRT(RL) 

DBPP(I)  =  DD( 11,2) 

2100  CONTINUE 
C 

I  =  IL- 1 

IU  =  IEND(l) 

DO  2200  11=2, IU 
I  =1  +  1 

DLS(I)  =  DELSTR( 11,1) 

VN( I )  =  DDD( II , 1)/SQRT(RL) 

DBPP(I)  =  DD( 11,1) 

2200  CONTINUE 

C  WRITE  (  6,200  )  (  DLS(I),  1=1, NT  ) 

C  WRITE  (  6,220)  (VN(I)  ,  1=1, NT) 

C 

RETURN 
END 

DATA  SET  KCBCCOMPGI  AT  LEVEL  001  AS  OF  08/24/84 

DATA  SET  KCBCCOMPGI  AT  LEVEL  001  AS  OF  08/24/84 

DATA  SET  KCBCCOMPGI  AT  LEVEL  003  AS  OF  08/24/84 

SUBROUTINE  COMPGI 

CALCULATE  GI 

COMMON  /BLCO/  NX,NXT,NP,NPT,NTR, IT, INVRS ,NS , IP 
COMMON  /BLC7/  C( 100,100) ,D( 100) ,DB( 100) ,DBP( 100) ,UE0( 100) ,GI 


180 

C 


260 

C 


SUM2 

= 

0.  0 

N1 

= 

NX  1 

N2 

= 

NXT 

DO  130  K 

= 

NT  ,N2 

SUM2 

= 

SUM2  + 

CONTINUE 

N1 

= 

2 

N2 

= 

NX-1 

SUM1 

= 

0. 

DO  260  K 

= 

N1,N2 

SUM1 

SUM1  + 

CONTINUE 

GI 

— 

UEO(NX) 

RETURN 

DBP(K) ) 


DBP(K) ) 

C(NX,NX)  *  DBP(NX) 


INT12240 
INT12250 
INT12260 
INT12270 
INT12280 
INT12290 
INT12300 
INT12310 
INT12320 
INT12330 
INT12340 
INT12350 
INT12360 
INT12370 
INT12380 
INT12390 
INT12400 
INT12410 
INT12420 
INT12430 
INT12440 
INT12450 
INT12460 
INT12470 
INT12480 
INT12490 
INT12500 
INT125 10 
INT12520 
INT12530 
INT12540 
INT12550 
INT12560 
INT12570 
INT12580 
INT12590 
INT12600 
INT12610 
INT12620 
INT12630 
INT12640 
INT12650 
INT12660 
INT12670 
INT12680 
INT12690 
INT12700 
INT127 10 
INT12720 
INTI 2 7 30 
INT12740 
INTI 2 750 
INT12760 
INT12770 
INT12780 


86 


i-^non  i~*ooo  o  o  o  nnoooooooooo 


r 

DATA  SET  KCBCDIFF3  AT  LEVEL  001  AS  OF  08/24/34 

DATA  SET  KCBCDIFF3  AT  LEVEL  001  AS  OF  08/24/84 

DATA  SET  KCBCDIFF3  AT  LEVEL  002  AS  OF  04/05/84 

SUBROUTINE  DIFF3  (N,X,F,FP,FPP,FPPP, IEND) 

DETERMINES  THE  DERIVATIVE  OF  THE  INPUT  FUNCTION  AT  THE  INPUT  PTS. 
DIMENSION  X( 101) ,F( 101) ,FP( 101) ,FPP( 101) ,FPPP( 101) 


FIRST  DERIVATIVES  USING  WEIGHTED  ANGLES 
N1=N- 1 

DX=Xf  2) -X( 1) 

DF=F( 2) -F( 1) 

ANG2=  ATAN2(DF,DX) 

DL2=DX 

DO  10  1=2, N1 
ANG1=ANG2 
DL1=DL2 
11=1+1 

DX=Xf Il)-X(I) 

DF=F( Il)-F(I) 

ANG2=  ATAN2(DF,DX) 

DL2=DX 

ANG=( DL2*ANG 1+DL 1*ANG2 ) / ( DL1+DL2 ) 

F?( I )=  TAN(ANG) 

IF  (I. NE. 2)  GO  TO  10 
ANGI  =  ANG 
ANG1I  =  ANGI 
CLI  =  DL1 

CONTINUE 

ANGF  =  ANG 

ANG2F  =  ANG 2 

DLF  =  DL2 

IEND1  =  IEND  +  1 

GO  TO  (11,12,13),  IEND1 

IEND  =  0  ,  EXTRAPOLATE  FOR  END  VALUES 

FP( 1 )  =  2. *(F( 2) -F( 1  ) )/DLI  -  FP(2) 

FP( N)  =  2.  *(Fi'N) -F(N1)  )/DLF  -  FP(N1) 

GO  TO  20 


IEND  =  1,  DERIVATIVES  ARE  CONTINUOUS  ACROSS  ENDS 

ANG  =  (DLI*ANG2F  +  DLF*ANG II)  /  (DLI  +  DLF) 

FP( 1 )  =  TAN ( ANG) 

FP(N)  =  FP( 1 ) 


INT12790 
INT12800 
INT12810 
INT12820 
INT12830 
INT12840 
INT12850 
INT12860 
INT12S70 
INT12880 
INT12890 
INT12900 
INT129 10 
INT12920 
INT12930 
INT12940 
INT12950 
INT12960 
INT12970 
INT12980 
INT12990 
INT13000 
INT13010 
INT13020 
INT13030 
INT13040 
INT13050 
INT13060 
INT13070 
INT13080 
INT13090 
INT13100 
INT13110 
INT13120 
INT13130 
INT13140 
INT13150 
INT13160 
INTI  3 170 
INT13180 
INT13190 
INT13200 
INT13210 
INT13220 
INT13230 
INT13240 
INT13250 
INT13260 
INT13270 
INT13280 
INT13290 
INT13300 
INT13310 
INT13320 
INT13330 


ST 


GO  TO  20 

n 

w 

C  I END  =  2,  IF  FIRST  DERIVATIVE  . LT.  0.0  RESET  TO  ZERO 
C 

13  CONTINUE 

FF( 1 )  =  2. *(F( 2) -F( 1  ))/DLI  -  FP(2) 

IF  (FP  (1)  .LT.  0.0)  FP  (1)  =  0.  0 
FP(N)  =  2.  *(F(N)-F(N1))/DLF  -  FP(N1) 

C 

C  SECOND  +  THIRD  DERIVATIVES  USING  CUBIC  FITS 

C 

20  DO  30  1=2, N1 

11  =1-1 

12  =1+1 

DX1  =  X  (II)  -  X  (I) 

DX2  =  X  (12)  -  X  (I) 

DF1  =2.0  *  ((F  (II)  -  F  (I))  /  DX1  -  FP  (I))  /  DX1 

DF2  =2.0  *  ((F  (12)  -  F  (I))  /  DX2  -  FP  (I))  /  DX2 

FPPP( I )=  3.0  *  (DF1  -  DF2)  /  (DX1  -  DX2) 

FP?  (I)=  DF1  -  DX1  *  FPPP  (I)  /  3.0 
30  CONTINUE 

FPPP( 1)=  FPPP  (2) 

FPPP(N)=  FPPP  (Nl) 

FPP  (1)=  FPP  (2  )  +  (X  (1)  -  X  (2  ))  *  FPPP  (2  ) 

FPP  (N)=  FPP  (Nl)  +  (X  (N)  -  X  (Nl))  *  FPPP  (Nl) 

C 

RETURN 

END 

C  DATA  SET  KCBCEDDY  AT  LEVEL  001  AS  OF  08/24/84 

C  DATA  SET  KCBCEDDY  AT  LEVEL  001  AS  OF  08/24/84 

C  DATA  SET  KCBCEDDY  AT  LEVEL  003  AS  OF  04/05/84 

SUBROUTINE  EDDY 
C 

C  CALCULATE  EDDY  VISCOSITY  USING  C  . S.  TWO-LAYERS  EDDY 
C  VISCOSITY  FORMULA 

COMMON  /BLCO/  NX,NXT,NP ,NPT,NTR, IT, INVRS ,NS , IP 
COMMON  /BLC1/  F( 101,2) ,U( 101 ,2) ,V( 101,2) ,W( 101,2) ,B( 101 ,2) 
C0MM0N/BLC7/  C( 100,100) ,D( 100) ,DB( 100) ,DBP( 100) ,UEO( 100) ,GI 
C0MM0N/EDDY1/  RL , RX , SQRX , RXNTR , GMTR , GMTRS ( 100) , 

+  ALFAS( 100) ,FFS( 100) ,RTS( 100) , IEDY,NXSPT 

COMMON  /GRD  /  ETA( 101) ,DETA( 101) , A( 101) 

COMMON  /GTY  /  X( 101) ,UE( 100) ,P1( 100) ,P2( 100) ,CEL,CELH 
DIMENSION  FINT(lOl) 

C 

C . 

c 

J0=1 

UDEL=0. 995*U(NP,2) 

DO  10  J=1 ,NP 
IF(U(J,2).  LT. UDEL)  JJ=J 
10  IF(U(J,2).  LT.  0.  0)  JO=J 

EDEL=ETA( JJ)+(ETA( JJ+1) -ETA( JJ) )/(U( JJ+1 , 2) -U( JJ, 2) ) 

+  *(UDEL-U( JJ,2) ) 

DO  15  J=1 ,NP 
ETADEL=ETA( J) /EDEL 
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IF( ETADEL.  GT.  1.  0)  ETADEL=1.  0 
15  FINT( J)=l.  0/(1.  0+5.  5 -"-ETADEL- *6) 

CALL  AMEAN( 1 , JJ ,ETA,FINT, 2) 

RU  =  RL 

IF  (IT  . GT.  1)  GO  TO  20 
ALFAS(N'X)  =  ALFAS(N'X-l) 

FFS(NX)  =  FFS(NX-l) 

RTS(NX)  =  RTS(NX-l) 

C 

GMTR  =  GMTRS(NX) 

IF  (NX  .  LE.  NS)  RU  =  RL  *  UE(NX) 

RL2  =  SQRT( RU  *  X(KX)) 

RL4  =  SQRT(RL2) 

RL216  =  0.  16  *  RL2 

20  VMAX  =  0.5  *  (V( 1,2)  +  V(l,l)) 

no  30  J=2 ,NP 

VM  =  0.5  *  (V(J,2)+V(J,1)) 

IF(VM  .GT.  VMAX)  VMAX  =  VM 
30  CONTINUE 

IF  ( IEDY  .EQ.  0)  GO  TO  35 

IF  (IT  .LE.  1  .OR.  GMTR  .  LT.  0.85  .OR.  NX  .  LE.  NTR+3  ) 

+  GO  TO  35 
C 

C  MODIFY  ALFA  USING  SIMPSON'S  ARGUMENTS 

C  CALL  SMPSON 

35  ALFA  =  ALFAS(NX) 

ED VO  =  ALFA  *  RL2  *  GMTR  *  (U(NP, 2)*ETA(NP)  -  F(NP,2)) 

DO  40  J=2 , NP 
JJ  =  J 

YBA  =  RL4*SQRT(VMAX)/26.  0*ETA(J) 

EL  =  1.  0 

IF(YBA  .LT.  10.0)  EL  =  1.  0  -  EXP(-YBA) 

EDVI  =  RL216  *  GMTR  *  ( EL*ETA( J) )**2  *  ABS(V(J,2)) 

IF( EDVI  .GT.  EDVO)  GO  TO  70 
B( J,2)  =  1.0  +  EDVI*FINT( J) 

IF(B( J,2)  .LT.  B( J-l , 2) )  B(J,2)  =  B(J-1,2) 

C  B(  J,  2)  =  1.0  +  EDVI 

40  CONTINUE 

JM  =2 

EM  =  B( 2 , 2) 

DO  50  J=2 ,NP 

IF(BM. GT.  B( J,2) )  GOTO  50 

JM  =  J 

BM  =  B(J,2) 

50  CONTINUE 

GOTO  90 

70  DO  80  J=JJ,NP 

80  B(J,2)  =  1.0  +  EDVO*FINT( J) 

C  80  B(J,2)  =1.0  +EDVO 
C 

90  CONTINUE 

BC  1,2)  =  1.0 

C 

JJ  =1 

DO  100  J=2 , NP 
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on  non  ooo  >-■  oooono  oi—  non 


IF(U(J,2)  .  LT.  0.  0)  JJ  =  J 
100  CONTINUE 

IF( JJ. EQ. 1)  GO  TO  110 

IN  THE  SEPARATED  REGION,  EDDY  VISCOSITY  IS  SET  EQUAL  TO 
THAT  IN  THE  PREV0US  STATION  TO  AVOID  NUMERICAL  TROUBLE 
JJP3  =  JJ  +  3 

JJP3  =  MIN( JJP3 ,  NP) 

CALL  AMEAN( 1 , J JP3 , ETA , B ( 1 , 2 ) , 2 ) 

10  CALL  AMEAN( 1 ,NP ,ETA, B( 1 , 2) , 1) 

RETURN 

END 

DATA  SET  KCBCEXTRAP  AT  LEVEL  001  AS  OF  08/24/84 

DATA  SET  KCBCEXTRAP  AT  LEVEL  001  AS  OF  08/24/84 

DATA  SET  KCBCEXTRAP  AT  LEVEL  008  AS  OF  02/13/84 

SUBROUTINE  EXTRAP ( NX, NXTE ,X,Y) 

EXTRAPOLATE  DATA  USING  PARABOLIC  FORMULA 
DIMENSION  X( 101) ,Y( 101) 


Y1  =  Y(NX-2) 

Y2  =  Y(NX-l) 

XI  =  X(NX-2) 

X2  =  X(NX-l) 

X3  =  X(NXTE) 

XI  =  XI  -X3 

X2  =  X2  -X3 

BB  =  (Y1-Y2)/(X1**2  -  X2**2) 

AA  =  Y1  -  BB  *  XI** 2 

DO  10  I=NX,NXTE 

XI  =  X(I)  -X3 

Y( I )  =  AA  +  BB  *  XI** 2 

CONTINUE 

RETURN 

END 

DATA  SET  KCBCFILLUP  AT  LEVEL  001  AS  OF  08/24/84 

DATA  SET  KCBCFILLUP  AT  LEVEL  001  AS  OF  08/24/84 

DATA  SET  KCBCFILLUP  AT  LEVEL  007  AS  OF  04/05/84 

SUBROUTINE  FILLUP( INDEX) 

COMMON  /BLCO/  NX , NXT , NP , NPT , NTR , IT , INVRS , NS , IP 

COMMON  /BLC1/  F( 101,2) ,U( 101,2) ,V( 101,2) ,W( 101,2) ,B( 101,2) 

COMMON  /GRD  /  ETAf 101) ,DETA( 101) ,A( 101) 


IF(NP.  GE.NPT  )  RETURN 
IF( INDEX. EQ. 2)  GOTO  10 

DEFINE  PROFILES  FOR  B.  L.  GROWTH 
NP1  =  NP  +  1 

NP  =  NP  +  2 

NP  =  MIN(NP,  NPT) 

NM  =  NP 

GOTO  20 
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C  FILL  UP  PROFILES  BEFORE  MOVING  TO  THE  NEXT  STATION 
C 


10 

NP1 

=  NP  +  1 

NM 

=  NPT 

20 

DO  30  : 

[=NP1  ,NM 

F( J . 2 ) 

=  F( J-l ,2)  +  DETA( J-1)*U( J-l ,2) 

U(  J,2) 

=  UCJ-1,2) 

V(J,2) 

=  0.  0 

B(-1, 2) 

=  B( j-  :  ,2) 

W  ( J ,  2  ) 

=  W( J-1,2) 

30  CONTINUE 


C 

RETURN 

END 

C  DATA  SET  KCBCHEADER  AT  LEVEL  001  AS  OF  08/24/84 

C  DATA  SET  KCBCHEADER  AT  LEVEL  001  AS  OF  08/24/84 

C  DATA  SET  KCBCHEADER  AT  LEVEL  001  AS  OF  04/05/84 

SUBROUTINE  HEADER  (  TITLE , SURFID , ISTRP  ) 

COMMON  /I SURF/  ISF 
C 

DIMENSION  TITLE(20),  SURFID(4) 

C 

10  FORMAT  (  1H1 ,2QX,20A4  ) 

20  FORMAT  (  1H0 , 15X, ’ BOUNDARY  LAYER  CALCULATION  FOR  '  , 

+  'UPPER  SURFACE  ' , 10X, ’ ICYCLE=’ , 15  /  16X,71(1H-)  /  ) 

30  FORMAT  (  1H0, 15X, ' BOUNDARY  LAYER  CALCULATION  FOR  ’ , 

+  'LOWER  SURFACE  ' , 10X , ' ICYCLE=' , 15  /  16X,71(1H-)  /  ) 

C 

C . 

c 

WRITE  (  6,10  )  TITLE 
IF( ISF  .EQ.  1)  WRITE  (  6,20  )  ISTRP 
IF( ISF  .EQ.  2)  WRITE  (  6,30  )  ISTRP 
C 

RETURN 

END 

C  DATA  SET  KCBCINPUT  AT  LEVEL  001  AS  OF  08/24/84 

C  DATA  SET  KCBCINPUT  AT  LEVEL  001  AS  OF  08/24/84 

C  DATA  SET  KCBCINPUT  AT  LEVEL  009  AS  OF  08/24/84 

SUBROUTINE  INPUT( ITR , ISWPT , SURFID) 

COMMON  /BLCO/  NX,NXT,NP,NPT,NTR,IT,INVRS,NS,IP 

COMMON  /BLC1 /  F( 101,2) ,U( 101,2) ,V( 101,2) ,W( 101 ,2)  ,B( 101,2) 

COMMON  /BLC2/  DELF( 101) ,DELU( 101) ,DELV( 101) ,DELW( 101) 

COMMON  /BLC7/  C( 100,100) ,D( 100) ,DB( 100) ,DBP( 100) ,UEO( 100) ,GI 
COMMON  /BONV/  ITMAX,EPSL,EPST,CONV 
C0MM0N/EDDY1/  RL,RX,SQRX,RXNTR,GMTR,GMTRS( 100) , 

+  ALFASf 100) ,FFS( 100) ,RTS( 100) , IEDY,NXSPT 

COMMON  /GTY  /  X( 101) ,UE( 100) ,P1( 100) ,P2( 100) ,CEL,CELH 
COMMON  /BLIN/  TITLE( 20) ,XC( 100) , YC( 100) , ISG( 100) ,DELS( 100) , 

+  XCTR , XTR , I STRP , ICYCLE , ICYTL , XCTRS( 2 ) ,TRFIND( 2) 

COMMON/TRN/  PGAMTR , OMEGA , RTHETB , RTRANB 
COMMON  /I SURF/  ISF 
DIMENSION  D 1 ( 100) ,D2( 100) ,D3( 100) 

DIMENSION  SUR1 ID(4) ,XCS( 100) 

LOGICAL  TRFIND 
C 
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c 

c 

INT15560 

INT15570 

ITMAX  =  15 

INT15580 

EPSL  =  0. 0001 

INT15590 

EPST  =  0. 01 

INT15600 

NPT  =  101 

INT15610 

ETAE  =  8.  0 

INT15620 

VGP  =  1.  14 

INT15630 

DETA1  =  0.  01 

INT15640 

NSS  =  NXT  /  4 

INT15650 

C  SEARCH  FOR  PRESSURE  PEAK  AS  SWITCH  POINT 

INT15660 

UMAX  =  UE( 1) 

INT15670 

DO  50  I  =  2  ,  NXT 

INT15680 

IF  (UE( I)  . LE.  UMAX)  GO  TO  55 

INT15690 

UMAX  =  UE( I) 

INT15700 

50 

CONTINUE 

INT157 10 

GO  TO  60 

INT15720 

55 

NS  =  I  -  4 

INT15730 

60 

IF  (NS  . GT.  NSS)  NS  =  NSS 

INT15740 

IF  (NS  .  LT.  3)  NS  =  3 

I.NT15750 

C 

INT157dO 

C 

CALCULATE  THE  PRESSURE  PARAMETERS  PI  +  P2  FOR  B.  L.  CALCULATION 

INT15770 

C 

USING  TRANSFORMED  COORDINATES 

INT15780 

C 

INT15790 

CALL  DIFF3  (NXT,  X,  UE,  Dl,  D2,  D3 .  0  ) 

INT15800 

DO  65  I  =  2, NXT 

INT15810 

P2( I)  =  X(I)  *  D1(I)  /UE( I) 

INT15820 

P1(I)  =  0.5  *  (1.0  +  P2( I ) ) 

INT15830 

65 

CONTINUE 

INT15840 

Pl(l)  =  0.5  *  (1.0  +  P2(l)) 

INT15850 

XCMIN  =XC(1) 

INT15860 

MIN  =  1 

INT15870 

DO  70  1=1, NXT 

INT15880 

IF( XCMIN.  LT.  XC( I) )  GOTO  70 

INT15890 

XCMIN  =  XC(I) 

INT15900 

MIN  =  I 

INT159 10 

70 

CONTINUE 

INT15920 

DO  80  I  =  1  ,  NXT 

INT15930 

IF  (I  .GE.  MIN)  THEN 

INT15940 

XCS(I)  =  XC(I) 

INT15950 

ISG(I)  =  1 

INT15960 

ELSE 

INT15970 

XCS(I)  =  -XCS(I) 

INT15980 

ISG(I)  =  -1 

INT15990 

END  IF 

INT16000 

80 

CONTINUE 

INT16010 

INVRS  =  NS  +  1 

INT16020 

C 

INT16030 

C 

SEARCH  FOR  TRANSITION  LOCATION 

INTI 6040 

ITRP1  =  ITR  +  1 

INT16050 

GOTO  (150,  95,  120,  150),  ITRP1 

INT16060 

C 

INT16070 

C 

TRANSITON  LOCATION  IS  TNPUT  (  =  XCTR) 

INT160S0 

95 

DO  100  1=1, NXT 

INT16090 

IF(XCTR.  LT.XCS(I))  GOTO  105 

INT16100 

100 

CONTINUE 

INT16110 

92 


NTH  =  N'XT  +  1 

INT16120 

GOTO  200 

INT16130 

105 

NTS  =1-1 

INTI 6 140 

IF  ( NTR.  LT.  3)  THEN 

INT16150 

NTR  =  3 

INT16160 

XCTR  =  XC(NTR) 

INT16170 

END  IF 

INT16180 

GOTO  200 

INTi 6190 

C 

INT16200 

C 

TRANSITION  LOCATION  IS  SET  AT  THE  PRESSURE  PEAK 

INT16210 

c 

INT16220 

120 

/'—N 

t-H 

V-' 

w 

'z> 

II 

s: 

INT16230 

IM  =  1 

INT16240 

DO  75  I  =  1 , NXT 

INT16250 

IF(UM. GT. UE( I) )  GOTO  75 

INT16260 

IM  =  I 

INT16270 

UM  =  UE(IM) 

INT16280 

75 

CONTINUE 

INT16290 

IF( IM. LT.  4)  IM  =  4 

INT16300 

NTR  =  IM 

INT16310 

XCTR  =  XC(NTR) 

INT16320 

GOTO  200 

INT16330 

C 

INT16340 

C 

TRAN ST I ON  LOCATION  WILL  BE  CALCULATED  BASED  ON  MICHEL  CRITERION 

INT16350 

c 

INT16360 

150 

NTR  =  NXT  +  1 

INT16370 

200 

DO  210  1=1, NXT 

INT16380 

210 

GMTRS( I )=  0.0 

INT16390 

C 

INT16400 

C  ' 

TRANSITION  LOCATION  PROVISIONALLY  FROM  PREVIOUS  CYCLE 

INT16410 

c 

INTI  6420 

IF  (TRFIND(ISl)  )  THEN 

INT16430 

DO  211  I  =  1  ,  NXT 

INT16440 

XCS(I)  =  XO' T) 

INT16450 

IF  (I  .Li.  MIN)  XCS(I)  =  -XCS(I) 

INT16460 

211 

CONTINUE 

INT16470 

DO  215  1=1, NXT 

INT16480 

IF(XCTRS(ISF)  .LE.XCS(I))  GOTO  217 

INT16490 

215 

CONTINUE 

INT16500 

217 

NTR  =  1-1 

INTI 65 10 

XCTR  =  XCTRS(ISF) 

INT16520 

IF  (NTR  .LT.  3)  THEN 

INT16530 

NTR  =  3 

INT16540 

XCTR  =  XC(NTR) 

INT16550 

END  IF 

INT16560 

END  IF 

INT16570 

C 

INT16580 

C  1 

CALCULATE  GAMTR  DISTRIBUTION 

INT16590 

c 

INT16600 

IF  (NTR.  GT.NXT-1)  GOTO  250 

INT16610 

FAC  =  (XCTR-XC(NTR) )/(XC(NTR+l) -XC(NTR) ) 

INT16620 

XTR  =  X(NTR)  +  FAC*(X(NTR+1) -X(NTR) ) 

INT16630 

UETR  =  UE(NTR)  +  FAC*(UE(NTR+1) -UE(NTR) ) 

INT16640 

RXNTR  =  XTR  *  UETR  *  RL 

INT16650 

GGFT  =  1 .  0 / PG AMTR* R L** 2 / RXNTR** 1.34 

INT16660 

DO  220  I=NTR+1,NXT 

INT16670 
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ALFAS(I)  =  0. 0168 
220  GMTRS(  I  )=  1.0 

ALFAS(NTR)  =  0.0168 
I'EINTG  =  0.0 
U1  =  0.  5/L'ETR 
XI  =  XTR 

DO  230  I=NTR+1 ,NXT 
U2  =  0.  5/UE( I) 

X2  =  X(I) 

UEINTG  =  UEINTG+(U1+U2)*(X2-X1) 

U1  =  U2 
XI  =  X2 

GG  =  GGFT*UEINTG*(X( I ) -XTR)*UE( I )**3 
IF(GG  . GT.  10. 0)  GO  TO  250 
GMTRS(I)  =  l.O-EXP(-GG) 

230  CONTINUE 
250  CONTINUE 
C 

C  GENERATE  B.  L.  ETA  GRIDS  +  INTIAL  VELOCITY  PROFILES 
C 

CALL  INTL( ETAE , DETA1 , VGP) 

DO  260  1=1, NXT 
UEO(I)  =  UE( I) 

260  CONTINUE 
C  PRINT  OUT  INPUT  DATA 
C 

IF  (ICYCLE  .GE.  ICYTL-1  .OR.  IP  .  GE.  0)  THEN 
CALL  HEADER(  TITLE , SURFID, ISTRP  ) 

WRITE( 6 , 1002)  NXT, ITR, IP, NS ,NTR, ISWPT 
VR ITE ( 6 , 100 3 )  VGP , DETA 1 , RL , XCTR , OMEGA , PGAMTR 
ELSE 

IF  (ISF.EQ. 1)  WRITE( 6 , 1008)  ICYCLE 
IF  (ISF.EQ. 2)  WRITE( 6 , 1009)  ICYCLE 
END  IF 

IF  (NTR.  LT.NXT)  THEN 

IF  (ITR.EQ. 1)  WRITE  (6,1005)  XCTR, XTR, NTR 
IF  (ITR. EQ.  2)  WRITE  (6,1006)  XCTR, XTR, NTR 
IF  (TRFIND(ISF))  WRITE(6, 1007)  XCTR, XTR, NTR 
END  IF 
RETURN 
C 

C  . . 

c 

2  F0RMAT(20A4) 

3  FORMATQOI5) 
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INTI 6 7 60 
INT16770 
INT16780 
INT16790 
INT16800 
INT168 10 
INT16820 
INT16830 
INT16840 
INT16850 
INT16S60 
INT16870 
INT16880 
INT16890 
INT16900 
INT169 10 
INT16920 
INT16930 
INT16940 
INT16950 
INT16960 
INT16970 
INT16980 
INT16990 
INT17000 
INTI 7 010 
INTI  7020 
INT17030 
INTI  7040 
INT17050 
INTI 7060 
INT17070 
INT17080 
INT17090 
INTI  7 100 
INT17110 
INT17120 


4  F0RMAT( 6F10. 0)  .  INT17130 

1001  F0RMAT( 1H1 , 20X, 20A4)  INT17140 

1002  FORMAT( 1H0 , 10H  NXT  =  ,I5,7X,10H  ITR  =  ,I5,7X/  INT17150 

+  1H  , 10H  IP  =  , 15 , 7X, 10H  NS  =  ,I5,7X/  INT17160 

+  1H  , 10H  NTR  =  , 15 , 7X, 10H  ISWPT=  ,15)  INT17170 

1003  FORMAT( 1H0 , 10H  VGP  =  ,E12.4,10H  DETA1=  ,E12.4/  INT17180 

+  1H  , 10H  RL  =  ,E12. 4 , 10H  XCTR  =  ,E12. 4/  INT17190 

+  1H  , 10H  OMEGA  =  ,E12.4,10H  PGAMTR=  ,E12.4)  INT17200 

1004  FORMAT( 1H0 , 3X , 2H  I , 6X , 2HXC , 1 IX , 2HYC , 1 IX , 2H  X , 1 IX , 2HUE , 1 IX , 2HP 1 ,  INT17210 

+  11X, 2HP2 , /( 1H  , 3X , 13 , 6E13. 5) )  INT17220 

1005  F0RMAT(/3X,’ BEGIN  OF  TRANSITION  IS  BEING  INPUT  AT  X/C  =',F8.4,4X,  INT17230 
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1006 


1007 


1008 

1009 

c 

c 

c 


c 

c 

c 

c 

c 


10 

20 


30 


40 

C 

C 

80 


90 


120 


+  'S/C  =' ,FS.4,4X,’NTR  =' ,13/) 

FORMAT(/3X, 'BEGIN  OF  TRANSITION  IS  SET  AT  PRESSURE  PEAK,  X/C  =  ’ , 
+  F?  4,4X, 1  S/C  =  '  ,F8.  4 ,4X , ' NTR  =',I3/) 

FORMAT^ /3X, ' BEGIN  OF  TRANSITION  IS  PROVISIONALLY  TAKEN  FROM 
+  'PREVIOUS  CYCLE:  X/C= ’ ,F8. 4 ,4X, ’ S/C  =’ ,F8. 4 ,4X, ’ NTR  =’ , 13/) 
FORMAT(//3X,’ UPPER  SURFACE  CALCULATIONS  OF  CYCLE1 ,13) 
FORMAT(//3X, ’LOWER  SURFACE  CALCULATIONS  OF  CYCLE', 13) 

END 

DATA  SET  KCBCINTL  AT  LEVEL  001  AS  OF  08/24/84 

DATA  SET  KCBCINTL  AT  LEVEL  001  AS  OF  08/24/84 

DATA  SET  KCBCINTL  AT  LEVEL  009  AS  OF  02/22/84 

SUBROUTINE  INTL(ETAE ,DETA1 ,VGP) 


INTI  7 240 
INTI 7 250 
INT17260 
INT17270 
INTI  7280 
INT17290 
INT17200 
INTI  73 10 
INT17320 
INT17330 
INT17340 
INT17350 


COMMON  /BLCO/  NX ,NXT,NP ,NPT,NTR , IT, INVRS ,NS , IP  INT17360 

COMMON  /ELC1/  F( 101 ,2) ,U( 101 ,2) ,V( 101 ,2) ,W( 101 ,2) ,B( 101 ,2)  INT17370 

COMMON  /BLC2/  DELF( 101) ,DELU( 101) ,DELV( 101) ,DELW( 101)  INT17380 

COMMON  /BLC6/  Sl( 101) ,S2( 101)  ,S3( 101)  ,S4( 101) ,S5( 101) ,S6( 101) ,  INT17390 

+  S7( 101), S8( 101), Rl( 101), R2( 101), R3( 101), R4( 101)1 NT 17400 

COMMON  /BCNV/  ITMAX,EPSL,EPST,CONV  INT17410 

COMMON  /GRD  /  ETA( 101) ,DETA( 101) ,A( 101)  INT17420 

COMMON  /GTY  /  X( 101) ,UE( 100) ,P1( 100) ,P2( 100) ,CEL,CELH  INT17430 

I NT 17 440 

.  I NT 17450 

INT17460 


GENERATE  THE  GRID 


INT17470 


INTI  7480 


DETA(l)  =  DETA1 

IF(  VGP.  LT.  1.  0)  VGP  =  1.  0 

IF(  ( VGP-1.  0)  .LE.  0.001)  GOTO  10 

NP  =  ALOG( (ETAE/DETA( 1))*( VGP-1. 0)+l. 0)/ALOG(VGP)  +  1.001 
GO  TO  20 

NP  =  ETAE/DETA(1)  +  1.001 
IFCNP  .LE.  NPT)  GO  TO  30 
WRITE( 6 ,  150  ) 

STOP 

ETA ( 1 )  =  0.  0 


DO  40  J=2 ,NPT 

ETA(J)  =  ETA(J-l)  +  DETA(J-l) 

DETAC  J)=  VGP"DETA( J-l) 

A(J)  =  0. 5*DETA(J-1) 

CONTINUE 

GENERATE  INITIAL  VELOCITY  PROFILE 
DO  90  J=1 ,NP 

ETAB  =  ETA( J)/ETA(NP) 

ETAB2  =  ETAB**2 

F(J,2)  =  0. 25*ETA(NP)*ETAB2*(3. 0  - 
U(J,2)  =  0. 5-ETAB*(3. 0  -  ETAB2) 

V ( J , 2 )  =  1.  5*( 1. 0  -  ETAB2)/ETA(NP) 

B( J,2)  =  1.  0 

W  ( J ,  2  )  =  1.  0 

CONTINUE 

NX  =1 

IT  =0 

IT  =  IT  +  1 

I F ( IT  .LT.  ITMAX)  GO  TO  130 


INT17490 
INTI  7500 
INTI 75 10 
INT17520 
INT17530 
INTI  7540 
INT17550 
INTI  7560 
INTI 75 70 
INT17580 
INT17590 
INT17600 
INTI  76 10 
INT17620 
INT17630 
INT17640 
INTI 7 650 
INTI  7660 
INT17670 
INTI  7680 

0. 5*ETAB2)  INT17690 

INTI  7 700 
INTI  7 7 10 
INT17720 
INT17730 
INTI 7 740 
INT17750 
INTI  7 7 60 
INT17770 
INTI 7 7 80 


oonnoooonoooo  ooo 


STOP 

130  CONTINUE 
C 

DO  140  J=  2, \P 

FB  =  0.5*(F(J,2)  +  F( J-l ,2) ) 

UB  =  0. 5*(U( J,2)  +  U( J-l ,2) ) 

VB  =  0. 5*( V( J , 2 )  +  V( J-1,2)) 

USB  =  0. 5*(U( J,2)**2  i  U( J-l , 2)**2) 

DERBV  =  (B(J,2)*V(J,2)  -  B( J- 1 ,2)*V( J-l ,2) )/DETA( J-l) 
FVB  =  0. 5*(F(J,2)*V(J,2)  +  F( J-l, 2)*V( J-1,2)) 

C 

S1(J)  =  0. 5*P1(NX)*F(J,2)  +  B( J,2)/DETA( J-l) 

S2( J)  =  0. 5*P1(NX)*F( J-1,2)  -  B( J-l ,2)/DETA( J-l) 

S3(J)  =  0. 5*P1(NX)*V(J,2) 

S4( J)  =  0.  5*P1(NX)*V(J-1,2) 

S5( J)  =  -P2(NX)*U( J, 2) 

S6( J)  =  -P2(NX)*U( J- 1,2) 

CRB  =  -P2(NX) 

R2( J)  =  CRB  -  (DERBV  +  P1(NX)*FVB  -  P2(NX)*USB) 

C 

Rl( J)  =  F( J-l ,2)  -  F(J,2)  +  DETA(J-1)*UB 
R3( J-l)=  U( J-1,2)  -  U( J,2)  +  DETA( J-1)*VB 
140  CONTINUE 

R 1  ( 1 )  =  0.  0 

R2( 1)  =  0. 0 

R3(NP)  =0.0 
CALL  S0LV3 

IF( ABS(DELV( 1) )  .  GT.  EPSL  )  GO  TO  120 
CALL  FILLUP( 2) 

CALL  OUTPUT(l) 

C 

RETURN 

C 

150  FORMAT( 1H0 , 37HNP  EXCEEDED  NPT  --  PROGRAM  TERMINATED) 

END 

DATA  SET  KCBCINTRP3  AT  LEVEL  001  AS  OF  08/24/84 

DATA  SET  KCBCINTRP3  AT  LEVEL  001  AS  OF  08/24/84 

DATA  SET  KCBCINTRP3  AT  LEVEL  003  AS  OF  04/05/84 

SUBROUTINE  INTRP3  (N1 ,X1 ,F1 ,FP1 ,FPP1 ,FPPP1 ,N2 ,X2 ,F2) 

CUBIC  INTERPOLATION 

GIVEN  THE  VALUES  OF  A  FUNCTION  (FI)  AND  ITS  DERIVATIVES 
AT  N1  VALUES  OF  THE  INDEPENDENT  VARIABLE  (XI) 

FIND  THE  VALUES  OF  THE  FUNCTION  (F2) 

AT  N2  VALUES  OF  THE  INDEPENDENT  VARIABLE  (X2) 

X2  CAN  BE  IN  ARBITRARY  ORDER 


DIMENSION  X1(101),F1(101),FP1(101) ,FPP1( 101) ,FPPP1( 101) ,X2( 101) , 
+  F2( 101) 

DATA  EPS  / IE-04/ 

C 


INT17790 
INI  17800 
INTI 78 10 
INTI  7820 
INT17830 
INTI  7840 
INT17850 
INT17S60 
INT17870 
INTI  7880 
INT17890 
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INT17920 
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INT17950 
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1 


JT  =  2 

DO  10  I  =  1  ,N'2 
XM  =  X2( I ) 

DO  20  J  =  JT,N'l 
J1  =  J  -1 
IF  (Xl(J).GE.  XM 
20  CONTINUE 

J  =  N1 
30  JT  =  J 

DXX 
DXX2 
DXX  3 

0  F2( I ) 


)  GO  TO  30 


=  X2( I )  -  X1(J1) 

=  DXX  *  DXX  /  2. 

=  DXX2  *  DXX  /  3. 

=  F1(J1)  +  DXX*FP1( Jl) 


+  DXX2*FPP1( Jl)  +  DXX3*FPPP1( Jl) 


RETURN 

END 

DATA  SET  KCBCJOIN  AT  LEVEL  001  AS  OF  08/24/34 

DATA  SET  KCBCJOIN  AT  LEVEL  001  AS  OF  08/24/84 

DATA  SET  KCBCJOIN  AT  LEVEL  012  AS  OF  02/20/84 

SUBROUTINE  J0IN( INDEX) 

COMMON  /BLCO/  NX , NXT , NP , NPT.NTR, IT, I NVRS , NS , I P 

F( 101,2) ,U( 101 ,2) ,V( 101 ,2) ,W( 101 ,2) ,B( 101 ,2) 

C( 100 , 100) , DC  100) ,DB( 100) ,DBP( 100) ,UE0( 100) ,GI 
RL , RX , SQRX , RXNTR , GMTR , GMTRS( 100 ) 

,ALFAS( 100) ,FFS( 100) ,RTS( 100) , IEDY.NXSPT 
X(101) ,UE( 100) ,P1( 100) ,P2( 100) ,CEL,CELH 
ETA( 101) ,DETA( 101) ,A( 101) 
FS(101),US(101),VS(101),BS(101),WS(101) 

VW( 100) , ITP( 100) , ISL( 100) ,DLS( 100) ,CF( 100) , 
THT( 100) ,NPSTR( 100) 


COMMON  /BLC1/ 
COMMON  /BLC7/ 
COMMON/EDDY 1/ 

h 

COMMON  /GTY  / 
COMMON  /GRD  / 
COMMON  /SAVE/ 
COMMON  /SMRY/ 


INDEX  =  1  FOR  THE  FIRST  SWEEP 
INDEX  =  2  FOR  SUBSEQUENT  SWEEP 

CALL  COMPGI 
CII  =  C(NX,NX) 

UES  *  GI  /  (1.0  •  DLS(NX)  *  SQRT(TL)  *  CII  ) 
IF( INDEX. EQ. 1)  GOTO  15 
C 

C  RETRIEVE  PROFILES  AT  STATION  NS  FOR  INVERSE  B.  L. 
C  CALCULATION 

DO  10  J=1 ,NPT 
F(J,2)  =  FS( J) 

U( J,2)  =  US(J) 

V( J, 2)  =  VS( J) 

W( J, 2)  =  WS(J) 

B( J,2)  =  BS( J) 

10  CONTINUE 

UES  =  UES/W( 1,2) 

SQS  =1.0 

GOTO  30 
15  CONTINUE 

SQS  =  1.0  /  SQRT(UES) 

DO  20  J=2 ,NPT 
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INT18370 

INT18380 

INT18390 

INT18400 

INT18410  ! 

INT18420 

INT18430 

INT18440 

INT18450 

INT18460 

INT18470 

INT16480 

INT18490 

INT18500 
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INT18600 
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INT18800 

INT18810 
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INT18860 

INT18870 
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INT1S890 

INT18900 
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20 

C 

30 


60 


C 


C 

C 


65 

70 


80 


C 

C 

C 

C 

C 


ETA(J)  =  ETA(.*1*SQS 
DETA(J-l)  =  E.A(J)  -  ETA(J-l) 

A(J)  =  0. 5*DETA(J-1) 

CONTINUE 

DO  60  J=1 ,NPT 
U  ( J ,  2 )  =  U(J,2)*UES 
W ( J , 2 )  =  UES  *  W ( J , 2 ) 

F(  J ,  2)  =  F  ( J ,  2  )  ■'■"  S  Q  S  ”  UES 
V(J,2)  =  V( J  j  2 ) / SQS*UES 
CONTINUE 
UE(NX)  =  W( 1 ,2) 

RX  =  UE ( NX ) *X( NX) *RL 

SQRX  =  SQRT(RX) 

IF(NX.  GT.  NTR)  CALL  EDDY 
CALL  FILLUP( 2) 

IF( INDEX. EQ. 2)  GOTO  70 

STORE  PROFILES  AT  STATION  NS  FOR  THE  NEXT  SWEEP 
DO  65  J=1 ,NPT 
FS(J')  =  F(  J,2) 

US(J)  =  U ( J , 2 ) 

VS( J)  =  V ( J , 2 ) 

WS(J)  =  W ( J , 2 ) 

BS(J)  =  B( J,2) 

CONTINUE 
DO  80  J=1 ,NPT 
F( J, 1 )  =  F( J , 2 ) 

U( J, 1)  =  U(J,2) 

V  ( J ,  1 )  =  V  ( J ,  2 ) 

W(J,1)  =  W(J,2) 

BC J, 1)  =  B( J,2) 

CONTINUE 

RETURN 

END 

DATA  SET  KCBCMAIN  AT  LEVEL  005  AS  OF  09/18/84 
PROGRAM  MAIN 


SUBROUTINE  CASBLP( K2 , XP , YP , XMP , YMP , XS , YS , DLSP , VC , DBPP 
+  ,RN,NBL, ITRI ,XCTRI ,CASEID) 

COMMON  /BLCO/  NX ,NXT,NP ,NPT,NTR , IT, INVRS ,NS , IP 
COMMON/BLIN/  TITLE ( 20) ,XC( 100) ,YC( 100) ,ISG( 100) ,DELS( 100) , 

+  XCTR , XTR , I STRP , ICYCLE , ICYTL , XCTRS (2) ,TRFIND(2) 

COMMON/EDDY 1/RL,RX, SQRX, RXNTR ,GMTR ,GMTRS( 100) , 

+  ALFAS( 100) ,FFS( 100) ,RTS( 100) , IEDY.NXSPT 

COMMON/BLOW/  VN(IOO) 

COMMON/TRN/  PGAMTR , OMEGA , RTHETB , RTRANB 
COMMON/PLOT/NVP( 2) ,NXVP( 20 , 2) , ICC 


INTI 89 10 
INT18920 
INT18930 
INT18940 
INT18950 
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INT18970 
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INT18990 
INT19000 
INT19010 
INT19020 
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INT19080 
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INT19120 
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INT19150 
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INT19180 
INT19190 
INT19200 
I  y'1'  19210 
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INT19240 
INT19250 
INT19260 
INT19270 
INT19280 
INT19290 
INT19300 
INT19310 
INT19320 
INT19330 
INT19340 
INT19350 
INT19360 
INT19370 
INT19380 
INT19390 
INT19400 


DIMENSION  CASEID(  20  ),  XCTRI  (  2  ),  ITRI  (  2  )  INT19410 
DIMENSION  XP  (  100),  YP  (  100),  XMP  (  100)  INT19420 
DIMENSION  YMP  (  100),  VC  (  100),  SM  (  100)  INT19430 
DIMENSION  XS  (  100),  YS  (  100),  NBL  (  2  )  INT19440 
DIMENSION  VT  (  100),  S  (  100),  DLSP  (  100)  INT19450 
DIMENSION  DLS  (  100),  XO  (  100),  YO  (  100)  INT19460 


DIMENSION  D1  (  100),  D2  (  100'* ,  D3  (  100) 

DIMENSION  XPS  (100),  YPS(IOO)  ,  DBPP(IOO) 

C 

10  FORMAT  (  20A4  ) 

20  FORMAT  (  415  ) 

30  FORMAT  (  3I5.3F10.0  ) 

40  FORMAT  (  6F10. 5  ) 

100  FORMAT  (  1H1 , 5X , 20A4  ) 

110  FORMAT  (  1H0,6X, ’CYCLE  NO.  =  ’,15  ) 

120  FORMAT  (  1H0 , 6X, ’ FINISHED  TOTAL  NUMBER  OF  CYCLES  =  ’,15, 

+  4X  ’ JOB  COMPLETED. ’  ) 

130  FORMAT  (  1H0 , f-X,  THE  UPDATED  DISPLACEMENT  SURFACE’ , /IX, 2HNX, 

+  9X . 2KXP , 12X , 2HYP , 1 IX , 3HDLS , 10X , 4HDBPP , 12X, 2HVN) 

140  FORMAT  (  I5,5E14. 6) 

150  FORMAT  (  1H0.6X,’READ  IN  CONTROL  POINTS  DATA 1 , / IX , 2HNX , 9X , 

+  3HXMP , 1 IX , 3HYMP , 1 2X , 2HSM , 1 3X , 2HVC ) 

160  FORMAT  (  I5.4E14. 6) 

170  FORMAT  (  1H0.6X, ’ INTERPOLATED  ORIGINAL  GEOMETRY  DATA’ , /IX, 

+  2HNX , 9X , 2HXP , 1 2X , 2HYP , 12X , 1HS , 14X , 2HVT) 

ISO  FORMAT  (I5.4E14. 6) 

C 

C . 

C 

GRANG(X1,X2,X3,Y1,Y2,Y3,X0)=  (X0-X2)*(X0-X3)/(X1-X2)/(X1-X3)*Y1 
+  +(X0-X1)*(X0-X3)/(X2-X1)/(X2-X3)*Y2+(X0-X1)*(X0-X2) 

+  /(X3-X1)/(X3-X2)*Y3 

ISWPT  =  1 
IEDY  =  1 
NBL(l)  =  91 

NBL( 2)  =  71 

C 

C  WRITE  (  6,100  )  CASEID 

5  CONTINUE 

ISTRP  =  ICYCLE 
NN  =  K2  -  1 
NT  =  K2 

C  WRITE  (  6,110  )  ICYCLE 

C  INTERPOLATE  OUTPUT  CONTROL  POINTS  TO  ORIGINAL  GEOMETRY 
C 

DO  15  I  *  1  ,  NT 
XPS(I)  =  XP(I) 

YPS(I)  =  YP(I) 

15  CONTINUE 

S(  1)  =  0.  0 

DO  50  I  =  2  ,  NT 

SCI)  =  s(l-l)  +  SQRT((XP(I)-XP(I-1))**2  + 

+  (YP(I)-YP(I-1))**2) 

50  CONTINUE 

SM( 1 )  =  SQRT((XMP(1)-XP(1))**2  +  ( YMP( 1) -YP( 1) )**2) 

DO  60  I  =  2  ,  NN 

SM(  I )  =  SM(I-l)  +  SQRT((XMP(  I)  -XMP(  I  - 1)  ),v*2+ 

+  ( YMP( I ) -YMP( I -1) )**2) 

60  CONTINUE 

C  CALL  AMEAN( NN- 10 , NN , SM , VC , 1 ) 

CALL  AMEAN( 1 , NN , SM , VC , 1 ) 

SNT  =  S( NT) 
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SM(NT)  =  SM(NN)+SQRT(  (XMP(NN)  -XP(N'T)  )**2+(  YMP(NN)  -YP(NT)  )**2) 
SMNT  =  S(  NT)  /  SM(N'T) 

DO  65  I  =  1  ,  NN 
65  SM(  I )  =  SM(  I )  SMNT 

CALL  DIFF3(NN, SM ,VC ,D1 ,D2 ,D3 , 0) 

CALL  INTRF3 ( NN , SM , VC ,D1 ,D2 ,D3 , NT , S , VT) 

C  PRINT  OCT  INPUT  DATA 
C 

C  WRITE (6 ,150) 

C  WRITE( 6 , 160)  ( I ,XMP( I) ,YMP( I) , SM( I) ,  VC( I) , 1=1 ,NN) 

C  WRITE( 6,170) 

C  WRITE( 6 , 180)  ( I ,XP( I),YP(I),S(I) ,  VT( I) , 1=1 , NT) 

C 

XPMIN  =  XP(1) 

DO  44  I  =  2  ,  NT 

IF  (XPC I)  .GT.  XPMIN)  GO  TO  44 

XPMIN  =  XP(I) 

44  CONTINUE 

CHORD  =  XP(NT)  -  XPMIN 

DO  45  I  =  1  ,  NT 

XP(I)  =  (XP( I ) -XPMIN)  /  CHORD 

YP( I )  =  YP(I)  /  CHORD 

45  CONTINUE 

CALL  COMPBL  (  CASEID ,XP , YP , VT, S ,DLSP ,DLS ,DBPP ,NBL, ITRI ,XCTRI , 
+  RN,NT, ISWPT) 

C 

C 

C  SMOOTH  THE  CALCULATED  DISPLACEMENT  THINKNESS 
C 

C  CALL  SMFIT(1,NT,S,DLS,D1,2) 

C 

C  ADD  DISPLACEMENT  THINKNESS  ON  THE  ORIGINAL  BODY 
C 

DO  70  I  =  1  ,  NT 
DLSP(I)  =  DLS(I) 

70  CONTINUE 

C  CALL  SMFIT  ( 1 ,NT,XS , YS ,D1 , 2) 

DO  80  I  =  1  ,  NT 
XP(I)  =  XPS(I) 

YP( I )  =  YPS(I) 

80  CONTINUE 

IF  (ICYCLE  .GE.  ICYTL-1  .OR.  IP  .  GE.  0)  THEN 
C  WRITE  (6,130) 

C  WRITE  (6,140)  ( I ,XP( I ) , YP( I ) ,DLS( I) ,DBPP( I) , VN( I) , 1=1 , NT) 

WRITE  (  6,120  )  ICYCLE 
END  IF 
C 

RETURN 

END 

C  DATA  SET  KCBCMAIN2  AT  LEVEL  001  AS  OF  08/24/84 

C  DATA  SET  ..JBCMAIN2  AT  LEVEL  001  AS  OF  08/24/84 

C  DATA  SET  KCBCMATN2  AT  LEVEL  010  AS  OF  04/06/84 

SUBROUTINE  MAIN2( ITR, ISWPT, SURFID) 

COMMON  /BLCO/  NX , NXT , NP , NPT , NTR , IT , INVRS , NS , I P 

COMMON  /BLC1/  F( 101 ,2) ,U( 101,2) ,V( 101 ,2) ,W( 101,2) ,B( 101,2) 

COMMON  /BLC2/  DELF( 101) ,DELU( 101) ,DELV( 101) ,DELW( 101) 
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INT20250 

INT20260 
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INT20290 

INT20300 

INT20310 
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INT20330 

INT20340 

INT20350 

INT20360 

INT20370 
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c 

c 


10 


15 

20 

25 


30 


70 


C 


COMMON  /BLC 7/  C(  100 , 100'S  ,D(  100)  ,DB(  100)  ,DBP(  100)  ,UEO(  100)  ,GI 

COMMON  / BONV/  ITMAX , EPSL , EPST . CONY 

COMMON/ EDDY 1 /  RL , RX , SQRX , RXSTR , GMTR , GMTRS ( 100) 

+  , ALFAS(  100)  ,FFS(  100)  ,RTS(  100) , IEDY,NXSPT 

COMMON  /GTY  /  X( 101) ,UE( 100) ,P1( 100) ,P2( 100) .CEL.CELH 
COMMON  / SMRY /  YW( 100) ,ITP( 100) ,ISL( 100) ,DLS( 100) ,CF(100) , 

+  THT( 100) ,NPSTR( 100) 

COMMON  /BLIN/  TITLE(20) ,XC( 100) ,YC( 100) ,ISG( 100) ,DELS( 100) , 

+  XCTR,XTR, ISTRP , ICYCLE , ICYTL,XCTRS( 2) ,TRFIND(2) 

C0MM0N/BLC9/  UEB(IOO)  ,  CFS(IOO) 

COMMON/TRN/  PGAMTR , OMEGA , RTHETB , RTRANB 
COMMON  /I SURF/  ISF 
COMMON/ PLOT/ NVP( 2) ,NXVP( 20,2), ICC 
DIMENSION  SURFID(4) ,RTSS( 11) 

LOGICAL  SMOOTH  ,  SEPART  ,  HOMOPY 

DATA  RTSS/O.  0,0.  1 , 0.  2 , 0.  3 , 0.  4 , 0.  5 , 0.  6 , 0.  7 , 0.  8 ,0.  9 , 1.  0/ 


GRANG(X1,X2,X3,Y1,Y2,Y3,X0)=  (X0-X2)*(X0-X3)/(X1-X2)/(X1-X3)*Y1 
+  +( XO  -X 1  )••  ( XO  -X3  )  /  ( X2  -X 1 )  /  ( X2  -X3 )  *  Y2+(  XO  -X 1  )*(  XO  -X2 ) 

+  /(X3-X1)/(X3-X2)*Y3 

ISWP  =  0 

INDEX  =  1 

IGROWT  =  2 

NXSPT  =  NXT  +  1 
CALL  JOIN( INDEX) 

NXSTOP  =  NXT- 1 

IF  (NS  .GE.  NTR  )  GOTO  15 

ISWP  =  ISWP  +  1 

NX  =  NX  +  1 

KOMOPY  =  .FALSE. 

CEL  =  0.  5*(X(NX)+X(NX-1))/(X(NX)-X(NX-1)) 

P 1  ( NX )  =  0.  5 
P2(NX)  =  0. 0 
CELH  =  0. 5*CEL 
IT  =0 
CALL  COMPGI 
IGR0W=1 

IT  =  IT  +  1 

RX  =  UE ( NX ) *X ( NX )  *RL 

SQRX  =  SQRT(RX) 

IF( IT  .LE.  ITMAX)  GO  TO  80 
1F( HOMOPY)  GO  TO  72 
IRC  =  1 

RT  =  RTSS(IRC) 

HOMOPY  =  .  TRUE. 

UEREF  =  UEO(NX-l) 

UESAVE  =  UEO(NX) 

UEO(NX)  =  RT*UESAVE+(  1.  0-RT)*UEREF 
DO  61  J=1 ,NP 
F(J,2)  =  F( J, 1) 

U( J ,2)  =  U(J,1) 

V (  J ,  2 )  =  V(  J  ,  1 ) 

W ( J , 2 )  =  W(J,1) 
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B( J , 2 )  =  B( J, 1) 

61  CONTINUE 

GO  TO  30 
C 

72  N’XSTOP  =  NX  -  1 

CALL  AMEANC NS , NXSTOP , X , CF , 1 ) 

C  CALL  AMEANC  NS , NXSTOP , X ,  VW , 1 ) 

C  CALL  HEADERC  TITLE , SURF ID , ISTRP  ) 

WRITE  (6,  250  )  ISWP 

WRITE  (6,  260  )  (M,XC(M) ,X(M) ,CF(M) ,DLS(M) ,THT(M) ,UE(M) , 

+  (JEO(M) ,D(M) ,DB(M) ,GMTRS(M) , ITP(M) ,NPSTR(M) ,M=1 .NXSTOP) 
WRITE(6,  270  )  NX 
STOP 

80  CONTINUE 

IF (NX  .GT.  NTR)  GOTO  100 
C 

C  LAMINAR  FLOW  CALCULATION 
C 

CALL  COEF( GAMMA 1 , GAMMA2 ) 

CALL  S0LV4( GAMMA 1 , GAMMA2 ) 

UE(NX)  =  U(NP, 2) 

IF(ABS(DELV(1))  .  GT.  EPSL)  GO  TO  70 
C 

C  CHECK  ON  LAMINAR  FLOW  SEPARATION.  IF  SEPARATION  OCCURS,  ASSIGN  BEGIN 
C  OF  TRANSITION  TO  THAT  POINT  AND  RECOMPUTE  THE  CURRENT  STATION  NX 
C 

IF(V( 1,2). GT.  0.  0  .OR.  ITR. NE.  3)  GOTO  110 
CALL  TRNS(ICODE) 

GOTO  25 
C 

C  TURBULENT  FLOW  CALCULATION 
C 

100  CONTINUE 
CALL  EDDY 

CALL  C OE F ( GAMMA 1 , GAMMA 2 ) 

CALL  S0LV4( GAMMA 1 , GAMMA2 ) 

UE(NX)  =  U(NP,2) 

VM  =  AMAX1(V(1,2),1.  0) 

IF(ABS(DELV( 1)/VM)  .  GT.  EPST)  GO  TO  70 
110  CONTINUE 
C 

C  CHECK  FOR  B.  L.  GROWTH 
C 

IF(NP  .GE.  NP7)  GO  TO  120 

IF( ABS(V(NP, 2) )  .LT.  0.0005  .AND.  ABS( 1.  0-U(NP-2,2)/U(NP,2)) 

+  .  LT.  0. 0035.  OR. IGROW.  GT.  IGROWT)  GOTO  120 

CALL  FILLUP'' n 
IGROW=IGROW+l 
IT  -  =  1 

GO  TO  70 
C 

120  CONTINUE 

CALL  FILLUPC 2) 

CALL  OUTPUTC  2 ) 

IF(NX.  GE.NTR  .OR.  ITR.  EQ.  0)  GOTO  150 
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c 

c 

c 


c 

c 

c 

c 


150 


154 

c 

c 

c 

c 

c 

c 

160 


165 

C 

C 


800 

810 

170 

C 


IF( NX.  LT.  3  .OR.  ITR.  NE.  3)  GOTO  150 

CALCULATE  TRANSITION  LOCATION  USING  MICHEL  METHOD 

CALL  TRNS(ICODE) 

IF( I COLE. EQ. 0)  GOTO  150 

TRANSITION  OCCURS  BASED  ON  MICHEL  CRITERIOR  AT  STATION  NX 
RECALCULATE  B.  L.  AT  NX  STATION  ASSUMING  THE  FLOW  IS  TRANSITIONAL 

IT  =0 

I GROW  =  1 

GOTO  70 
CONTINUE 

IF( . NOT.  HOMOPY  )  GO  TO  154 
IF(  RT  .GT.  0.9999)  GO  TO  154 
IRC  =  IRC  +  1 
RT  =  RTSS(IRC) 

l:ec(nx’i  =rt*uesave  +  (i.  o-rt)*ueref 
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GO  TO  30 
CONTINUE 

IF(NX  .LT.  NXSTOP)  GO  TO  20 

THE  B.  L.  CALCULATION  FOR  THE  CURRENT  SWEEP  IS  COMPLETED. 
CHECK  FOR  THE  CONVERGENCE  AND  ,  IT  NOT,  MOVE  TO  THE  NEXT 
SWEEP. 


CONTINUE 

D(NXT)  =  GRANG(X(NXT-3) ,X(NXT-2) ,X(NXT-1) ,D(NXT-3) ,D(NXT-2)  , 

+  D(NXT-l) ,X(NXT)) 

DLS(NXT)=  GRANG( X( NXT- 3 ) ,X(NXT-2) ,X(NXT-1) ,DLS(NXT-3) ,DLS(NXT-2) , 
+  DLS(NXT-l) ,X(NXT) ) 

UE(NXT)  =  GRANG(X(NXT-3) ,X(NXT-2) ,X(NXT-1) ,UE(NXT-3) , 

+  UE(NXT-2),UE(NXT-1),X(NXT)) 

DO  165  I  =  1  ,  NXSTOP 
CFS(I)  =  CF( I ) 

CALL  AMEAN( NS, NXSTOP, X,CF,1) 

CALL  AMEANC NS , NXSTOP . X , VW , 1 ) 

CALL  HEADER*  TITLE , SURFID , ISTRP  ) 

IF( ICYCLE  .LT.  ICYTL-1  .AND.  IP  . LT.  0)GO  TO  170 
WRITE  (6,  250  )  ISWP 

WRITE  (6,  262  )  1 ,XC( 1) ,X( 1) ,CF( 1) ,DLS( 1) ,THT( 1) ,UE( 1) , 

+  UE0( 1) ,0.  0, GMTRS( 1 ) , ITP( 1 ) , NPSTR( 1 ) 

WRITE  (6,  264  )  (M,XC*M) ,X(M) ,CF(M) ,DLS(M) ,THT(M) ,UE(M) , 

+  UEO(M) ,DLS(M)/THT(M) ,GMTRS(M) , ITP(M) ,NPSTR(M) ,M=2, NXSTOP) 

IF  ( ( ICYCLE.  EQ. ICYTL). AND. ( IP. EQ. -2). AND. (NVP( ISF). NE. 0) )  THEN 
WRITE* 12,800)  NS+1 ,NTR,XCTR 

WRITER  12,810)  (XC(M) ,X(M) ,UE(M) ,CF(M) ,GMTRS(M) , ISG(M) , 

2  M=l, NXSTOP) 

FORMAT*. 215, F10.  6) 

FORMAT* 5E15.  5,15) 

END  IF 
CONTINUE 

DMAX  =  D( 1) 
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180 

C 

c 
c 
c 


DDMAX  =  ABS(D( 1)  -  DB(1)) 

DO  ISO  I  =  2 , NXT 

DMAX  =  AMAX1C  DMAX,D(I)  ) 

DD  =  ABS( D( I )  -  D3( I ) ) 

DDMAX  =  AMAX1C  DDMAX ,DD  ) 

CONTINUE 

IF  (  ABS(  DDMAX  /  DMAX  )  .  LE. 


UPDATE  D  FOR  THE  NEXT  SWEEP 


0. 0050  )  RETURN 


190 

195 

200 

205 


210 


250 


260 


262 


264 

270 


C 

C 

C 


IF  (ISWP  .GT.  1)  GO  TO  195 
DO  190  I  =  NS  ,  NXT 

D(  I)  =  D(  I)*(  1.  0-HOMEGA*(UE(  I)/UEO(  I ) -1.  0)  ) 

GO  TO  205 

IF  (ISWP  .EQ.  2)  GOTO  205 
DO  200  I  =  NS  ,  NXT 

D( I)  =  D( I)  *  (1. 0+OMEGA*(UE(I)/UEB(I)-1.  0)) 

IFCISWP  .GE.  ISWPT  )  RETURN 

NX  =  NS 

NP  =  NPSTRCNX1 

INDEX  =  2 

DO  210  1=  1 , NXT 

DB( I)  =  D( I ) 

UEB(I)  =  UE( I ) 

CONTINUE 
GOTO  10 

FORMATC 1H0 , '  **  SUMMARY  OF  INVERSE  BOUNDARY  LAYER  SOLUTIONS.  **' , 
+  1H0,4X,'ISWP  =' , 13/ ) 

FORMATC 1H0 , 4X , 2HNX , 5X , 3HX/C , 9X , 1HX , 9X , 2HCF ,  8X , 

+  3HDLS , 8X , 3HTHT , 9X , 2HUE ,  8X , 3HUE0 ,  10X , 1HD  ,9X,2HDB,3X, 

+  4HGMTR,4X,2HIT, lX,2HNP/( 1H  , 3X, 13 ,F10.  5 ,8E11.  4,F8.  4, 213) ) 

FORMATC 1H0 , 4X , 2HNX , 6X , 3HX/C , 1 IX , 1HX , 10X , 2HCF , 9X , 

+  3HDLS , 9X , 3HTHT , 10X , 2HUE , 9X , 3HUE 0 , 1 IX , 1HH , 8X , 

+  4HGMTR,4X,2HIT,1X,2HNP/C1H  , 3X, 13, 9E12. 4,213)) 

FORMATC 1H  , 3X, 13, 9E12.  4,213) 

FORMATC 1H0,'  **  ITERATIONS  EXCEEDED  ITMAX  AT  NX  =’ ,15, ’ ,**' ,/ 

+  1H0 , '  **  CALCULATIONS  STOP.  **') 

END 

DATA  SET  KCBCOUTPUT  AT  LEVEL  001  AS  OF  08/24/84 

DATA  SET  KCBCOUTPUT  AT  LEVEL  001  AS  OF  08/24/84 

DATA  SET  KCBCOUTPUT  AT  LEVEL  002  AS  OF  02/22/84 

SUBROUTINE  OUTPUTC INDEX) 

COMMON  /BLCO/  NX,NXT,NP,NPT,NTR,IT,INVRS,NS,IP 
COMMON/BLIN/  TITLEC 20) ,XCC 100) , YC( 100) , ISG( 100) ,DELS( 100) , 

+  XCTR , XTR , I STRP , I CYCLE , I C YTL , XCTRS  C  2 ) , TRF I ND ( 2 ) 

COMMON  /BLC1/  F( 101 ,2) ,UC 101,2) , VC  101,2) ,WC 101 ,2) ,B( 101,2) 

COMMON  /BLC7/  C( 100 , 100) ,DC 100) ,DB( 100) ,DBPC 100) ,UEO( 100) ,GI 
COMMON/EDDY 1 /  RL , RX , SQRX , RXNTR , GMTR , GMTRS C 100) 

+  ,ALFAS( 100) ,FFS( 100) ,RTSC 100) , IEDY,NXSPT 

COMMON  /GTY  /  X( 101) ,UE( 100) ,P1( 100) ,P2C 100) ,CEL,CELH 
COMMON  /GRD  /  ETAC 101) ,DETA( 101) ,AC 101) 

COMMON  /SMRY/  VWC 100) , ITPC 100) , ISLC 100) ,DLSC 100) ,CFC 100) ,THTC 100) 
+  NPSTR(IOO) 

COMMON  /I SURF/  ISF 
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COMMON/ PLOT/NVP( 2)  ,NXVP( 20 , 2) , ICC 
C 

C  . 

c 

IIP (NX)  =  IT 
NPSTR(NX)=NP 
IF(NX.  GT.  1)  GOTO  5 
DLS(NX)=  0.0 
VW(NX)  =  0.  0 
D(NX)  =0.0 
THT(NX)=  0.  0 
CF(NX)  =  0.  0 
YW( NX)  =  0.0 
GOTO  150 

5  GOTO  (10,100,200),  INDEX 
C 

C  CALCULATE  B.  L.  PARAMETERS  FOR  TRANSFORMED  COORDINATES 
10  CONTINUE 

CF(NX)  =  2.0  *  V( 1 , 2)  *  B( 1 , 2) /SQRX 

VW  ( NX )  =  UE(NX)*SQRT(UE(NX)/X(NX))*V(1,2) 

DLS(NX)=  X(NX) /SQRX  *  (ETA(NP)-F(NP,2)) 

D(NX)  =  UE(NX)  *  DLS(NX)  *  SQRT(RL) 

U1  =  U( 1 , 2)  *  (1. 0  -U( 1,2)) 

SUM  =  0.  0 
DO  20  J=2 ,NP 

U2  =  U(J,2)  *  (1.  0  -U( J, 2) ) 

SUM  =  SUM  +  A(J)  *  (U1  +  U2) 

U1  =  U2 

20  CONTINUE 

THT(NX)=  X( NX) /SQRX  *  SUM 
GOTO  150 
C 

C  CALCULATE  B.  L.  PARAMETERS  FOR  SEMI-TRANSF  COORDINATES 
100  CONTINUE 

SQXC  =  SQRT(X(NX) ) 

SQRL  =  SQRT(RL) 

CF(NX)  =  2. 0  *  V( 1 , 2)  *  B ( 1 , 2 ) / ( SQXC*SQRL*W ( NP , 2 ) **2 ) 
VW ( NX )  =  V(l,2)  /  SQXC 
UE(NX)  =  U(NP,2) 

RX  =  RL  *  UE(NX)  *  X(NX) 

DLS(NX)  =  (ETA(NP) -F(NP , 2)/U(NP, 2) )/SQRL*SQXC 
SUM  =0.0 

U1  =  U(1,2)/U(NP,2)*(1. 0  -U( 1 , 2)/U(NP , 2) ) 

DO  120  J=2 , NP 

U2  =  U(J,2)/U(NP,2)*(1.  0  -U(J,2)/U(NP,2)) 

SUM  =  SUM  +  A(J)  *  (U1  +  U2) 

U1  =  U2 

120  CONTINUE 

THT(NX)  =  SUM  /SQRL  *  SQXC 
D(NX)  =  (U(NP , 2)*ETA(NP) -F(NP, 2) )  *  SQXC 
150  IF  (NX  .GE.  NXT)  GO  TO  160 

IF  ( IEDY  .EQ.  0  .OR.  NX  .  LE.  NTR+2)  GO  TO  160 
C 

C  MODIFY  ALFA  USING  SIMPSON'S  ARGUMENTS 
C 
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160 


175 

C 

C 

200 


C 

210 


4001 

4000 

4100 

4200 

4300 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 


CALL  SMPSON 
DO  175  J=1 ,NPT 
F(J,1)  =  F(J,2) 

U(J,1)  =  U(J  2) 

Y(J,1)  =  V( J ,2) 

W(J,1)  =  V(J,2) 

B( J, 1)  =  B(J.2) 

CONTINUE 

IF  ((IP.  LE.  0).  AND.  ( (IP.  NE. -2).  OR.  (ICYCLE.LT.  ICYTL)))  RETURN 

PRINT  OUT  VELOCITY  PROFILES 
IF  (NX.  EQ.  1)  GOTO  210 
IF  (NX. LE.NS)  THEN 

FAC1  =  SQRT(X(NX)/RL/UE(NX)) 

FAC2  =  1.0 
ELSE 

FAC1  =  SQRT(X(NX)/RL) 

FAC2  =  1.  0/UE(NX) 

END  IF 


INT23370 
INT23380 
INT23390 
INT23400 
INT23410 
INT23420 
INT23430 
INT23440 
INT23450 
INT23460 
INT23470 
INT23480 
INT23490 
INT23500 
INT235 10 
INT23520 
INT23530 
INT23540 
INT23550 


NPM1  =  NP  -1 
WRITE( 6 ,4001)  NX,X(NX) 

WRITE( 6,4000) 

WRITE( 6,4100)  ( J,ETA( J) ,F( J,2) ,U( J,2) ,V( J,2) ,W( J,2) ,B( J,2) , 

+  ETA( J)*FAC1 ,U( J, 2)*FAC2 ,J=1 ,NPM1 ,3) 

WRITE( 6 , 4100)  NP,ETA(NP) ,F(NP,2) ,U(NP,2) , V(NP, 2) ,W(NP, 2) ,B(NP , 2) , 
+  ETA(NP)*FAC1,U(NP,2)*FAC2 

IF  (IP. NE.-2)  RETURN 

IF  ( (NXVP( ICC,ISF).  NE.  NX).  OR.  ( ICC.  GT.  NVP(ISF)))  RETURN 

WRITE( 12,4200)  NP 

WRITE( 12 , 4300)  (ETA( J) , J=1 ,NP) 

VRITE( 12,4300)  (U( J,2) , J=1,NP) 

ICC  =  ICC+1 
RETURN 

FORMAT( / 1H0 , ' NX  =',I5,'  S/C  =  ' ,F10. 5) 

FORMAT! 1H0 , 2H  J , 9X , 3HETA , 15X , 1HF , 13X , 1HU , 13X , 1HV , 13X , 1HW , 1 3X , 1HB , 
+  13X,3HY/C, 10X,4HU/UE) 

FORMAT! 1H  ,I3,E14.  5,2X,5E14.  5,2X,2E14. 5) 

FORMAT( 15) 

FORMAT! 8F 10. 6) 

END 

DATA  SET  KCBCSMFIT  AT  LEVEL  001  AS  OF  08/24/84 

DATA  SET  KCBCSMFIT  AT  LEVEL  001  AS  OF  08/24/84 

DATA  SET  KCBCSMFIT  AT  LEVEL  001  AS  OF  08/15/83 

SUBROUTINE  SMFIT(NS ,ND,X,Q,D,KS) 

THIS  SUBROUTINE  SMOOTHES  DATA,  Q,  USING  FIVE-POINT  FORMULA. 

NS  :  BEGINNING  POINT?  ND  :  END  POINT 
X  :  INDEDENPENT  COORDINATE?  D  :  WORKING  STORAGE 
Q  :  VARIABLE  TO  BE  SMOOTHED 
KS  :  NO  OF  SMOOTHING 

DIMENSION  X( 101) ,Q( 101) ,D( 101) 

SMT5(Q1 ,Q2 ,Q3 ,Q4 ,Q5  )  =  0. 0625*( 10. 0*Q3+4. 0*(Q2+Q4)-Q1-Q5) 

SMT3(Q1 ,Q2 ,Q3 ,X1 ,X2 ,X3)  =  0.  5*(Q2+(Q1*ABS(X3-X2)+Q3*ABS(X2-X1) ) 
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C 

C 

C 

20 

40 

100 

C 

200 

220 

250 

300 


+  /ABS(X3-X1) ) 

IF(KS.LE.  0)  RETURN 
NS PI  =  NS+1 

NSP2  =  NS+2 

NDM1  =  ND-1 

NDM2  =  ND-2 

NDIF  =  ND-NS+1 

IF  (  NDIF  .LT.  3  )  RETURN 

IF  (  NDIF  .  LT.  5  )  GO  TO  200 

DO  100  K=1 ,KS 

D(NS+1)=  SMT3fQ(NS) ,Q( NS+1) ,Q( NS+2) ,X(NS) ,X(NS+1) ,X(NS+2)) 
D(ND-1)=  SMT3-  XND-2) ,Q(ND-1) ,Q(ND) ,X(ND-2) ,X(ND-1) ,X(ND) ) 
DO  20  I=NSP2 ,NDM2 

D( I)  =  SMT5( Q( 1-2) ,Q( 1-1) ,Q( I) ,Q( 1+1) ,Q( 1+2) ) 

CONTINUE 

DO  40  I=NS:i ,NDM1 
Q(  I )  =  D( I) 

CONTINUE 

RETURN 

DO  300  K  =  1,KS 
DO  220  I  =  NSP1.NDM1 

D( I )  =  SMT3(Q( 1-1) ,Q( I) ,Q( 1+1) ,X( 1-1) ,X( I) ,X( 1+1) ) 

CONTINUE 

DO  250  I  =  NSP1 ,NDM1 
Q( I)  =  D(I) 

CONTINUE 

CONTINUE 

RETURN 

END 

DATA  SET  KCBCSOLV3  AT  LEVEL  001  AS  OF  08/24/84 

DATA  SET  KCBCSOLV3  AT  LEVEL  001  AS  OF  08/24/84 

DATA  SET  KCBCSOLV3  AT  LEVEL  005  AS  OF  02/21/84 

SUBROUTINE  SOLV3 


INT23930 

INT23940 

INT23950 

INT23960 

INT23970 

INT23980 

INT23990 

INT24000 

INT24010 

INT24020 

INT24030 

INT24040 

INT24050 

INT24060 

INT24070 

INT24080 

INT24090 

INT24100 

INT24110' 

INT24120 

INT24130 

INT24140 

INT24150 

INT24160 

INT24170 

INT24180 

INT24190 

INT24200 

INT24210 

INT24220 

INT24230 

INT24240 

INT242S0 

INT24260 

INT24270 

INT24280 

INT24290 

INT24300 


COMMON  /BLCO/  NX , NXT , NP , NPT , NTR , IT , INVRS , NS , IP  INT24310 

COMMON  /BLC1/  F( 101 ,2) ,U( 101,2) ,V( 101,2) ,W( 101,2) ,B( 101,2)  INT24320 

COMMON  /BLC2/  DELF( 101) ,DELU( 101) ,DELV( 101) ,DELW( 101)  INT24330 

COMMON  /BLC6/  Sl( 101) ,S2( 101) ,S3( 101) ,S4( 101) ,S5( 101) ,S6( 101) ,  INT2>340 

+  S7(101),S8(101) ,R1( 101) ,R2C 101) ,R3( 101) ,R4( 101) INT24350 

COMMON  /GRD  /  ETA( 101) ,DETA( 101) ,A( 101)  INT24360 

COMMON  /BLCB/  All( 101) ,A12( 101) ,A13( 101) ,A14( 101) ,  INT24370 

+  A21( 101) ,A22( 101) ,A23( 101) ,A24( 101)  INT24380 

.  INT24390 


A 1 1  ( 1)=  1.0 
A12( 1)=  0. 0 
A13(  1)=  0.0 
A21(  1)=  0.0 
A22(  1)=  1.0 
A23(  1)=  0.0 
Gll=-1.  0 


INT24400 

INT24410 

INT24420 

INT24430 

INT24440 

INT24450 

INT24460 


G12=-A( 2)  INT24470 

G13=  0.0  INT24480 
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G2 1=  S4( 2) 

G23=-S2( 2) / A( 2) 

G22=  G23+S6C  2 ) 

FORWARD  SWEEP 

DO  500  J=2,NP 

IF(J  . EQ.  2)  GO  TO  100 

DEN'  =  (A13( J-1)*A21( J-l) -A23( J-1)*A11( J-l) -A( J)* 

+  (A12( J-1)*A21( J-l)-A22( J-1)*A11( J-l) ) ) 

DENI  =  A22( J-1)*A( J) -A23( J-l) 

Gll=  ( A23( J-l)+A( J)*(A( J)*A21( J-l) -A22( J-l) ) ) /DEN 
G12=  -( A(J)*A(J)+G11*(A12(J-1)*A(J)-A13( J-l) )) /DENI 
G13=  (G11*A13(J-1)+G12*A23(J-1))/A(J) 

G21=  (S2( J)*A21( J-l) -S4( J)*A23( J-l)+A( J)*(S4( J)* 

+  A22( J-l) -S6( J)*A21( J-l) ) )/DEN 

G22=  ( -S2( J)+S6( J)*A( J) -G21*(A(J)*A12( J-l) -A13( J-l) )) /DENI 
G23=  G21*A12( J-1)+G22*A22( J-l) -S6( J) 

3  All(  J)=  1.  0 

A12(J)=-A(J)-G13 
A13( J)=  A(J)*G13 
A21( J)=  S3( J) 

A22( J)=  S5( J) -G23 
A23( J)=  S1(J)+A(J)*G23 

Rl( J)  =  R1(J)-(G11*R1(J-1)+G12*R2(J-1)+G13*R3(J-1)) 

R2(J)  =  R2( J)-rG21*Rl( J-1)+G22*R2( J-1)+G23*R3( J-l) ) 

3  CONTINUE 


BACKWARD  SWEEP 

DELUi'NP)  =  F 
El  =  F 
E2  =  F 
DELV(NP)  =  ( 


R3(NP) 

R1(NP)-A12(NP)*DELU(NP) 

R2(NP) -A22(NP)*DELU(NP) 

(E2*A11(NP)-E1*A21(NP))/(A23(NP)*A11(NP)-A13(NP)* 
A21(NP) ) 

(E1-A13(NP)*DELV(NP) )/All(NP) 


DELF(NP)  *  (E1-A13(NP)*DELV(NP) )/All(NP) 

J  =  NP 

J  =  J-l 

E3  =  R3( J)-DELU( J+l)+A( J+1)*DELV( J+l) 

DEN 2  =  A21(J)*A12(J)*A(J+1)-A21(J)*A13(J)-A(J+1)*A22(J)*A11(J)+ 

(-  A23( J)*A11( J) 

DELV(J)  =  (A11(J)*(R2(J)+E3*A22(J))-A21(J)*R1(J)-E3*A21(J)*A12(J) 
y  )/DEN2 

DELU(J)  =-A( J+1)*DELV( J) -E3 

DELF(J)  =  (Rl( J)-A12( J)*DELU( J)-A13( J)*DELV( J) )/All( J) 

IF(  J  .GT.  1)  GO  TO  600 


DO  700  J=1,NP 

F( J, 2)=  F( J,2)+DELF( J) 

U(J,2)=  U(J,2)+DELU(J) 

V( J, 2)=  V(J,2)+DELV(J) 

CONTINUE 
U(  1 , 2)=  0.0 

CALL  EDGCHK( NP ,ETA , F( 1 , 2) , U( 1 , 2) , V( 1 , 2) ) 
RETURN 


INT24500 
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INT25000 
INT25010 
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DATA  SET  KC8CS0LV4  AT  LEVEL  001  AS  OF  08/24/84 

DATA  SET  KCBCS0LV4  AT  LEVEL  001  AS  OF  08/24/84 

DATA  SET  KCECS0LV4  AT  LEVEL  001  AS  OF  02/21/84 

SUBROUTINE  S0LV4( GAMMA 1 , GAMMA2 ) 

COMMON  /ELC0/  NX , NXT , NP , NPT , NTR , IT , INVRS , NS , I P 

COMMON  /BLC1/  F( 101 ,2) ,U( 101,2) ,V( 101 ,2) ,W( 101,2) ,B( 101 ,2) 

COMMON  /BLC2/  DELF( 101) ,DELU( 101) ,DELV( 101) ,DELW( 101) 

COMMON  /BLC6/  Sl( 101)  ,S2( 101) ,S3( 101) ,S4( 101) ,S5( 101) ,S6( 101) , 
l-  S7(  101)  ,S8(  101)  ,R1(  101)  ,R2(  101)  ,R3(  101)  ,R4( 

COMMON  /GRD  /  ETA( 101) ,DETA( 101) ,A( 101) 

COMMON  /BLCB/  All( 101) ,A12( 101) ,A13( 101) ,A14( 101) , 
i-  A21(  101)  ,  A22(  101)  ,  A23(  101)  ,A24(101) 


A 1 1 C 1 )  =  1.  0 

A12(  1 )  =  0.  0 

A 13 ( 1 )  =  0.  0 

A14( 1)  =  0. 0 

A21( 1)  =  0. 0 

A22(l)  =  1.  0 

A23(  1)  =  0. 0 

A24(  1)  =  0. 0 

DO  10  J  =  2,NP 

AA1  =  A13(J-1)-A(J)*A12(J-1) 

AA2  =  A23( J-l) -A( J)*A22( J-l) 

AA3  =  S2(J)-A(J)*S6(J) 

DET  =  AA2*A11(J-1)-AA1*A21(J-1) 

AJS  =  A( J)**2 

Gil  =  -(AA2+A21f  J-1)*AJS)/DET 

G12  =  (A11(J-1)*AJS+AA1)/DET 

G13  =  A12(J-1)*G11+A22(J-1)*G12+A(J) 

G14  =  A14(J-1)*G11+A24(J-1)*G12 

G21  =  ( S4( J)*AA2 -A2 1 ( J* 1 )*AA3 ) /DET 

G22  =  (A11(.T-1)*AA3-S4(  J)*AA1)/DET 

C-23  =  A12(J  *)*G21+A22(J-1)*G22-S6(J) 

G24  =  A14( J- 1 )"G21+A24( J- 1 )*G22-S8( J) 

Ail(J)  =  1.  0 
A12(  J)  =  -G13 

A13( J)  =-  A(J)*G13 

A14( J)  =  -G14 

A21( J)  =  S3( J) 

A22( J)  =  S5( J) -G23 

A23( J)  =  Sl( J)+A( J)*G23 

A24( J)  =  S7( J)-G24 

Rl( J)  =  Rl( J)  -G11*R1(J-1)-G12*R2(J-1)-R3(J-1)*G13 

t-  -G14*R4(  J-l) 

R2(J)  =  R2( J)  -G21*R1( J-l) -G22*R2( J-l) -R3( J-1)*G23 

E  -G24*R4( J-l) 

CONTINUE 


BACKU'ARD  SWEEP 
J  =  NP 

G1  =  GAMMA 1/GAMMA2 

R3( J)  =  R3( J) /GAMMA2 
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Rl(  J)  =  R1(J)-A12(J)*(R4(J)+R3(J))-A14(J)*R3(J)  INT25610 

R2( J)  =  R2(J)-A22(J)*(R4(J)+R3(J))-A24(J)*R3(J)  INT2S620 

Cl  =  A11(J)-G1*(A12(J)+A14(J))  INT25630 

C2  =  A21(J)-G1*(A22(J)+A24(J))  INT25640 

DET  =  C1*A23(J)-C2*A13(J)  INT25650 

DELF(J)  =  (R1(J)*A23(J)-R2(J)*A13(J))/DET  INT25660 

DELV(J)  =  ( C1*R2( J) -C2*R1( J) )/DET  INT25670 

DELW(J)  =  R3( J) -G1*DELF( J)  INT25660 

DELU(J)  =  R4(J)+DELV(J)  INT25690 

20  J  =  J-l  INT25700 

CC1  =  DELU( J+l) -R3( J) -A( J+1)*DELV( J+l)  INT25710 

CC2  =  DELW( J+l ) -R4( J)  INT25720 

CC3  =  A13(J)-A(J+1)*A12(J)  INT25730 

CC4  =  Rl( J) -A12( J)*CC1-A14( J)*CC2  INT25740 

CC5  =  A23(J)-A(J+1)*A22(J)  IKT25750 

CC6  =  R2( J) -A22( J)*CC1-A24( J)*CC2  INT25760 

DENO  =  A11(J)*CC5-A21(J)*CC3  INT25770 

DELF(J)  =  ( CC4*CC5 -CC3*CC6 ) /DENO  INT25780 

DELV(J)  =  (All(  J)  •VCC6-A21(  J)*CC4)/DEN0  INT25790 

DELW(J)  =  CC2  INT25800 

DELU(J)  =  CC1-A(J+1)*DELV(J)  INT25810 

IF(  J  .  C-E.  2)  GO  TO  20  INT25820 

DO  30  J  =  1  ,NP  INT25830 

F( J,2)  =  F( J, 2)+DELF( J)  INT25840 

U( J,2)  =  U(J,2)+DELU(J)  INT25850 

V( J,2)  =  V(J,2)+DELV(J)  INT25860 

W(J,2)  =  W(J,2)+DELW(J)  INT25870 

30  CONTINUE  INT25880 

U( 1 , 2)  =  0.0  INT25890 

CALL  EDGCHK(NP,ETA,F(1,2),U(1,2),V(1,2))  INT25900 

RETURN  INT25910 

C  .  INT25920 

END  INT25930 

DATA  SET  KCBCTRNS  AT  LEVEL  001  AS  OF  08/24/84  INT25940 

DATA  SET  KCBCTRNS  AT  LEVEL  001  AS  OF  08/24/84  INT25950 

DATA  SET  KCBCTRNS  AT  LEVEL  005  AS  OF  03/13/84  INT25960 

SUBROUTINE  TRNS(ICODE)  INT25970 

INT25980 

CALCULATE  TRANSITION  LOCATION  USING  MICHEL  CRITERION  INT25990 

INT26000 

COMMON  /BLCO/  NX ,NXT,NP ,NPT,NTR, IT, INVRS ,NS , IP  INT26010 

COMMON  /BLC1/  F( 101 ,2) ,U( 101 ,2) ,V( 101,2) ,W( 101 ,2) ,B( 101,2)  INT26020 

COMMON/BLIN/  TITLE(20) ,XC( 100) ,YC( 100) ,ISG( 100) ,DELS( 100) ,  INT26030 

+  XCTR,XTR , ISTRP, ICYCLE , ICYTL,XCTRS( 2) ,TRFIND( 2)  INT26040 

COMMON/EDDY1/  RL,RX,SQRX,RXNTR,GMTR,GMTRS( 100)  INT26050 

+  , ALFAS( 100) ,FFS( 100) ,RTS( 100) ,IEDY,NXSPT  INT26060 

COMMON  /GRD  /  ETA( 101) ,DETA( 101) ,A( 101)  INT26070 

COMMON  /GTY  /  X( 101) ,UE( 100) ,P1( 100) ,P2( 100) ,CEL,CELH  INT26080 

COMMON /TRN/  PGAMTR , OMEGA , RTHETB , RTRANB  INT26090 

C  .  INT26100 

100  FORMAT(/3X. 'ELGIN  OF  TRANSITION  HAS  BEEN  DETECTED  BY  MICHEL"S  INT26110 

+  ' CRITERION:  X/C  =' ,F8. 4,4X, ' S/C  =' ,F8.  4,4X, 'NTR  =’,I3/)  INT26120 

110  FORMAT(/3X, 'BEGIN  OF  TRANSITION  IS  ASSUMED  AT  THE  POINT  OF  ' ,  INT26130 

+  ’ LAMINAR  SEPARATION:  X/C  =’ ,F8.  4,4X, ' S/C  =' ,F8. 4,4X, ' NTR  =' ,13/)  INT26140 

C  .  INT26150 

ICODE  =  0  INT26160 
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ISEP  =  1 
C 

IF(V(  1,2).  LT.  0.  0)  THEN’ 

C  ***  TRANSITION  PROCESS  HAS  BEGUN  DUE  TO  LAMINAR  SEPARATION  *** 

FAC  =  V(  1 , 1 )  /(  V(  1 , 1 )  -V(  1 , 2) ) 

GOTO  20 
END  IF 
C 

C  ***  CHECK  MICHEL'S  TRANSITION  CRITERION  *** 

ISEP  =  0 

SUM  =0.0 

FI  =  U( 1 , 2)/U(NP,2)*( 1.  0-U( 1 ,2)/U(NP,2) ) 

DO  10  J=2 ,NP 

F2  =  U(J,2)/U(NP,2)*(1. 0-U(J,2)/U(NP,2)) 

SUM  =  SUM  +  (FI  +  F2  )  *A(J) 

FI  =  F2 

10  CONTINUE 

CONV  =  SQRT( RL/X( NX) ) 

IF(NX. LE.NS)  CONV  =  SQRT(RX)/X(NX) 

THETA  =  SUM  /  CONV 

RTHETA  =  RL  *  UE(NX)  *  THETA 

RTRAN  =  1.174  *  ( 1. 0+22400. 0/RX)  *  RX**0.46 

IF(RTHETA.LT.  RTRAN)  THEN 
RTHETB  =  RTHETA 
RTRANB  =  RTRAN 
RETURN 
END  IF 
C 

C  ***  TRANSITION  PROCESS  HAS  BEGUN  BECAUSE  OF  MICHEL'S  CRITERION  *** 
FAC  =  ( RTHETB -RTRANB ) / ( RTRAN-RTRANB -RTHETA+RTHETB ) 

***  COMPUTE  EXACT  LOCATION  OF  TRANSITION  BEGIN  *** 

20  NTR  =  NX-1 
NTR 1  =  NTR  +  1 

XCTR  =  XC(NX-l)  +  FAC*(XC(NX)-XC(NX-1)) 

XTR  =  X(NX-l)  +  FAC"(X(NX) -X(NX-l) ) 

UETR  =  UE(NX-l)  +  FAC*(UE(NX) -UE(NX-l) ) 

IF  (  ISEP  .EQ.  0  )  WRITE  (6,100)  XCTR , XTR , NTR 
IF  (  ISEP  .EQ.  1  )  WRITE  (6,110)  XCTR, XTR, NTR 
ICODE  =  1 

***  CALCULATE  INTERMITTENCY  DISTRIBUTION  *** 

RXNTR  =  XTR  *  UETR  *  RL 
GGFT  =  RL** 2 / PGAMTR / RXNTR** 1 .  34*UETR**3 
DO  30  I=NTR1,NXT 
ALFAS(I)  =  0. 0168 
GMTRS(  I  )=  1.0 
30  CONTINUE 

ALFAS(NTR)  =  0.0168 
UEINTG  =0.0 
U1  =  0.  5/UETR 
XI  =  XTR 

DO  40  I=NTR1 ,NXT 
U2  =  0.5/UECI) 

X2  =  X( I) 

UEINTG  =  UEINTG+(U1+U2)*( X2-X1) 
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U1  =  U2 

INT26730 

XI  =  X2 

INT26740 

GG  =  GGFT*UE I NTG* ( X(  I )  -XTR) 

INT26750 

IF( GG  . GT.  10.0)  GOTO  50 

INT26760 

GMTRS(I)  =  l.O-EXP(-GG) 

INT26770 

40  CONTINUE 

INT267SO 

c 

INT26790 

c  - 

RESET  FINITE  DIFFERENCE  CALCULATIONS  *** 

INT26800 

50  DO  60  J=1 ,NPT 

INT26810 

F( J , 2)  =  F(J,1) 

INT26820 

U(J,2)  =  U(J,1) 

INT26830 

V  ( J ,  2  )  =  V(  J ,  1 ) 

INT26840 

B( J,2)  =  B( J, 1 ) 

INT26850 

W( J, 2)  =  K( J, 1) 

INT26860 

60 

CONTINUE 

INT26870 

RETURN 

INT26880 

END 

INT26890 

C 

INT26900 

c 

INT269 10 

SUBROUTINE  EDGCHK( NP ,  ETA,  F,  U,  V) 

INT26920 

c 

INT26930 

DIMENSION  ETA(lOl),  F(101),  U(101),  V(101) 

INT26940 

n 

INT26950 

JS  =  NP  -  3 

INT26960 

NPM1  =  NP  -  1 

INT26970 

DO  10  J=JS,  NPM1 

INT26980 

JJ  =  J 

INT26990 

IF(U( J). GE. U(NP)  .OR.  V( J) .  LT.  0.  0)  GOTO  20 

INT27000 

10 

CONTINUE 

INT27010 

RETURN 

INT27020 

20 

JS  =  JJ  -  1 

INT27030 

IF( JS. GT. (NP-2) )  JS  =  NP-2 

INT27040 

CALL  AMEAN( JS ,  NP,  ETA,  U,  1) 

INT27050 

CALL  AMEANCJS,  NP,  ETA,  F,  1) 

INT27060 

DETAP  =  ETA(JS)  -ETA(JS-l) 

INT27070 

VJP  =  (U( JS) -U(JS-l) ) /DETAP 

INT27080 

DO  30  J=JS ,NPM1 

INT27090 

DETAM  =  ETA( J+1)-ETA( J) 

INT27100 

VJM  =  (U(J+1)-U(J)) /DETAM 

INT27U0 

V(J)  =  (VJM*DETAP  +  V JP*DETAM ) / ( DETAP+DETAM ) 

INT27120 

VJP  =  VJM 

INT27130 

DETAP  =  DETAM 

INT27140 

30 

CONTINUE 

INT27150 

V(NP)  =  -V(NP-l)  +  2.  0  *  VJP 

INT27160 

RETURN 

INT27170 

C 

**-.V-V*-.V*-V****************V--V***V.-Vf*****V.-**Vc************Vr**-V********** 

INT27180 

C 

NOTES:  (FOR  CHANGING  FROM  THE  ORIGINAL  PROGRAM) 

* 

INT27190 

C 

k 

INT27200 

C 

1.  'EDDY'  HAS  BEEN  MODIFIED  BY  ADDING  ' FINT' . 

* 

INT27210 

C 

2.  SUBROUTINE  'EDGCHK'  HAS  BEEN  ADDED. 

k 

INT27220 

C 

3.  GROWTH  LIMIT  HAS  BEEN  ADDED  FOR  2  IN  'MAIN2'. 

* 

INT27230 

C 

* 

INT27240 

C 

* 

INT27250 

C 

*  VfVrVf-.VVc  *  iskivitltkirk  kkkkkk  k  kkkkkk  k  kkkk  kkk  kkkkkkkkk  kkkkkkkkkkitki 

tkick  k  kkitk 

INT27260 

END 

INT27270 
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SUBROUTINE  SMPSON 

COMMON  /BLCO/  NX ,NXT,NP,NPT,NTR , IT, INVRS ,NS , IP 
COMMON  /ELC1/  F( 101,2) ,U(101,2),V( 101 ,2) ,W( 101,2) ,B( 101,2) 
COMMON  / 3 LC 7 /  Cl.  100,100) , DC  100)  ,DB(  100)  ,DBP(  100)  ,UE0(  100)  ,GI 
COMMON/EDDY 1 /  RL , RX , SQRX , RXNTR , GMTR , GMTRSC 100 ) 

+  , ALFAS( 100) ,FFS( 100) ,RTS( 100) , IEDY.NXSPT 

COMMON  /GTY  /  X( 101) ,UE( 100) ,P1( 100) ,P2( 100) ,CEL,CELH 
COMMON  /GRD  /  ETA( 101) ,DETA( 101) ,A( 101) 

DIMENSION  CRD( 12) ,RTD( 12) 

DATA  RTD/O. 00,0. 05,0. 12,0. 20,0. 30,0. 40,0. 50,0. 60,0. 70, 

+  0.80,0.90,1.00/ 

D.iTA  CRD/5.  00,4.  75,4.  35 , 3.  80 , 3.  25 , 2.  85 , 2.  58 , 2.  37 , 2.  25  , 

+  2.15,2.06,2.00/ 


STEP  1  CALCULATE  (DU/DX)/(DU/DY) 

IF( NX. LT. NXSPT)  GOTO  10 

IN  THE  SEPARATED  REGION,  ALFA  SET  TO  BE  CONSTANT 
ALFAS(NX)=  ALFASP 
RETURN 

C( . 

C 

C  10  CONTINUE 

C  IFCVC 1,2). GT.  0. 0)  GOTO  20 

C 

C  SEPARATION  OCCURS.  ALFA  SET  TO  BE  THE  PREVIOUS  ITERATED  VALUE 
C  ALFASP  =  ALFAS(NX) 

C  NXSPT  =  NX 

C  RETURN 

C( . 

C  MODIFY  OUTER  EDDY  BASED  ON  SIMPSON  SUGGESTION 
TM  =0.0 
JM  =  1 

DO  30  J=2 ,NP 

TS  =  (B(J,2)-1.  0)*  V(J,2) 

IF(TS.LT.  TM)  GOTO  30 
TM  =  TS 

JM  =  J 

30  CONTINUE 

VNXM  =  0.  5*(V(JM,2)+VCJM,1)) 

IF  (NX  .  LE.  NS)  GOTO  35 

DUDX  =  (U( JM, 2) -U( JM, 1) )  /  (X(NX) -X(NX-l) ) 

GO  TO  38 

35  DUDX  =  CEL*(U(JM,2)-U(JM,1))+P2(NX)*U(JM,2)+0.5*ETA(JM)* 

+  VNXM* ( P2 ( NX) - 1 . 0) 

38  RU  =  RL 

IF(NX. LE. NS)RU  =  RL  *  UEO(NX)  *  X(NX) 

RL2  =  SQRT(RU) 

RR  =  DUDX/VNXM/RL2 

C 

C  STEP  2  :  CALCULATE  (UU  -  VV)/UV 
VNXM  =  0.  5*(V( 1 ,2)+V( 1,1)) 
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INT27420 
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INT27470 
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c 


c 

c 

c 

c 

c 

c 

60 

c 

c 

70 


C( 

c 


10 

20 

30 


40 


RT  =  VNXM/TM 

PRINT' (3X, 215, 3F10.  3) ’ , NX , JM , VNXM , TM , RT 
IF  (RT  . LT.  0.  0)  RT  =  0.  0 
IF( RT.  GT.  1. 0)  GOTO  60 

CR  =  6.0  /( 1.  0  +  2.0  *  RT*(  2.  0  -RT) ) 

OR  =2.0 

DO  40  1=2,12 

IF(RT. LT. RTD(i))  GOTO  50 
40  CONTINUE 
GOTO  70 

50  CR  =  CF.E(  I-i)+( CRD( I) -CRD( 1-1) )*(RT-RTD( 

GOTO  70 

CR  =  (1.  0  +  RT  )  /RT 
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STEP  3  :  CALCULATE  FF 
FR  =  CR  *  RR 
IF(FR  .  GT.  0.  35)  FR  =  0.  35 
IF  (FR  . LT.  -0.  8)  FR  =  -0.  8 
FFS(NX)=  (FFS(NX)  +  (1.0  -FR))/  2.0 

RTS(NX)=  RT 

ALFAS(NX)=  0.  0168/FFS(NX)**2.  5 
RETURN 


END 

SUBROUTINE  XSPACE ( NI , NRITE , XI I ,XLLT,RAD ,NL1 ,NR1) 
DIMENSION  XII( 200) ,T( 200) 

DATA  PI/3. 14159265359879/ 

RAD  =  PI 

NLEFT=NI-1 -NRITE 
NR4=NRITE/2 

IF((NRITE/2*2)  .  NE.  NRITE)  NR4=(NRITE+l)/2 

NL1=NR4+1 

NL2=NR4+NLEFT 

NR1=NL2+1 

NR2=NI 

PI2=0. 5*PI 

RAD2=( PI -RAD) /2. 0+PI2 

RAD3=RAD2+RAD 

SRT  =RAD2/FL0AT(NR4) 

SRT2=SRT 

IF( (NRITE/2*2)  .  NE.  NRITE)  SRT2=RAD2/FLOAT( NR4- 1 ) 
SLT  =  RAD/FLOAT(NLEFT) 

DO  10  1=1 ,NR4 

XII(I)=0. 5*( 1.  0+COS(FLOAT(I-1)*SRT)) 

DO  20  I=NL1,NL2 

XII(I)=0.  5+XLLT*COS(FLOAT( I-NL1)*SLT+RAD2) 

DO  30  I=NR1 ,NR2 

XII(I)=0.  5*( 1.  0+C0S(FL0AT( I-NR1)*SRT2+RAD3) ) 
NA=(NI+l)/2 

IF( (NI/2*2)  .EQ.  NI)  NA=NI/2+l 
FN1=FL0AT(NA-1) 

FN2=FN1 

IF( (NI/2*2)  .EQ.  NI)  FN2=FLOAT( NA-2) 

DO  40  1=1, NA 
T(I)=FL0AT(NA-I)/FN1 
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CALL  AMEAN’C  1 .NA.T.XII , 1) 

XDIF  =  XII(l)  -  X I I ( 2 ) 

IF( XDIF  .  LT.  0.004)  THEN 
DO  45  1=2,5 

XII(I)  =  XII(I-n-XDIF*3.  0 
CONTINUE 

CALL  AMEAN( 2,NA,T,XII,10) 

END  IF 

DO  50  I=NA,NI 
XII(I)  =  XIICNI-I+l) 

RETURN 

END 

SUBROUTINE  TRGRID  (N1  ,  XO  ,  YO,NI  ,NRITE ,XLLT,N10 ,RAD , ID ,NXSS) 
THIS  SUB.  IS  TO  REGRID  SPACING  NEAR  TRAILING-EDGE 

DIMENSION  XO( 200) , YO( 200) ,XI( 200) ,YI( 200) ,D1( 200) ,D2( 200) ,D3( 200) 
+  X00( 200) , YOO( 200) ,XII(200) ,YII(200) ,WX(200) ,WY(200), 

+  WXI  200),WYI(200),T(200) 

N20  =  N10 

I F ( (N 1/2*2)  . FO.  Nl)  N20  =  N10+1 
IF(  (NI-(NI,'1)*2)  .NE.  0)  Nil®  (MI-l)/2+l 
IFUNI-(NI/2)*2)  .EQ.  0)  NlI=NI/2 
N2I  =  Nil 

IF C (NI/2*2)  .EQ.  NI)  N2I  =  N1I+1 

CALL  XSPACE(NI,NRITE,XI,XLLT,RAD,NL1,NR1) 

PRINT  * , ' NRITE= ' , NRITE , '  XLLT= ' , XLLT 
WRITE  (6,  290) 

WRITE  (6,  300)  (XO(I)  ,I=1,N1) 

WRITE  (6,  298) 

WRITE  (6,  300)  ( YOC I)  ,I=1,N1) 

IF( ID  .EQ.  2)  THEN 
DO  60  I=NL1 ,NXSS 
YI( I )=YO( I) 

XI(I)=XO(I) 

NXST=NXSS+8 

XM1=(XI(NXST) -XI(NXSS) )/8. 0 
XM2=(XI(NR1) -XI(NXSS) )/(NRl-NXSS) 

DO  62  I=.NXSS,NR1-1 

XI ( I )=(XO( I )+XI( I ) )/2. 0 

DO  65  I=NXSS ,NXST 

T( I )=FLOAT( I -NXSS)*XM1+XI(NXSS) 

CALL  AMEAN(NXSS,NXST,T,XI ,8) 

DO  68  I=NXSS ,NR1 
T(I)=FLOAT(I-NXSS)*XM2+XI(NXSS) 

CALL  AMEAN(NXSS ,NR1 ,T,XI ,12) 

N10=NL1 

N1I=N10 

N20=N1+1-NXSS 

N2I=NI-NXSS+1 

ELSE 

NXSS=N2I 

END  IF 

FOR  LOWER  SURFACE 
DO  5  I  =  1  ,  N10 
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II  =  N10  -  I  +  1 

INT28950 

VX(II)  =  XO(I) 

INT28960 

VY (II)  =  YO( I ) 

INT28970 

5 

CONTINUE 

INT28980 

DO  7  I  =  1  ,  Nil 

INT28990 

II  =  Nil  -  I  +  1 

INT29000 

WXI(II)  =  XI(I) 

INT29010 

7 

CONTINUE 

INT29020 

CALL  DIFF3(N10 ,WX,VY ,D1 ,D2 ,D3 , 0) 

INT29030 

CALL  INTRP3(N10,WX,WY,D1,D2,D3,N1I,WXI,WYI) 

INT29040 

DO  9  I  =  1  ,  Nil 

INT29050 

II  =  Nil  -  I  +  1 

INT29060 

Y I  ( II)  =  WYI(I) 

INT29070 

9 

CONTINUE 

INT29080 

C 

INT29090 

C 

FOR  UPPER  SURFACE 

INT29100 

DO  10  I  =  1  ,  N20 

INT29110 

II  =  N1  -  N20  +  I 

INT29120 

XOO(I)  =  XO(II) 

INT29130 

YOO(I)  =  YO(II) 

INT29140 

10 

CONTINUE 

INT29150 

DO  20  I  =  1  ,  N2I 

INT29160 

II  =  N1  -  N2I  +  I 

INT29170 

WXI(I)  =  XI(II) 

INT29180 

20 

CONTINUE 

INT29190 

CALL  DIFF3(N20,X00 , YOO ,D1 ,D2 ,D3 ,0) 

INT29200 

CALL  INTRP3 ( N20 , XOO , YOO , D 1 , D2 , D3 , N2 I , WXI , WYI ) 

INT29210 

DO  25  I  =  1  ,  N2I 

INT29220 

II  =  N2I  -  I  1 

INT29230 

XII(II)  =  VXI(I) 

INT29240 

YII(II)  =  WYI(I) 

INT29250 

25 

CONTINUE 

INT29260 

C 

INT29270 

C 

COMBINE  TWO  SURFACES  INTO  ONE  CIRCLE 

INT29280 

C 

INT29290 

C 

NN  =  N1  -  N10  -  N20 

INT29300 

C 

DO  30  I  =  1  ,  NN 

INT29310 

C 

II  =  Nil  +  I 

INT29320 

C 

I 11=  N10  +  I 

INT29330 

C 

XI(II)  =  XO(III) 

INT29340 

C 

YI(II)  =  YO(III) 

INT29350 

C30  CONTINUE 

INT29360 

DO  40  I  =  1  ,  N2I 

INT29370 

11  =  NXSS  +1-1 

INT29380 

12  =  N2I  -  I  +  1 

INT29390 

XI(I1)  =  XI 1(12) 

INT29400 

Y I  ( 11)  =  YII( 12) 

INT29410 

40 

CONTINUE 

INT29420 

XI(1)  =  X0(1) 

INT29430 

XI(I1)=  X0(N1) 

INT29440 

YI( 1)  =  Y0( 1) 

INT29450 

YI( Il)=  Y0(N1) 

INT29460 

C 

INT29470 

N1  =  11 

INT29480 

DO  50  I  =  1  ,  N1 

INT29490 

XO(I)  =  XI(I) 

INT29500 
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Y0(I)  =  YI( I) 

INT295 10 

50 

CONTINUE 

INT29520 

C 

INT29530 

c 

WRITE  (6  ,  295) 

INT29540 

c 

WRITE  (6  ,  300)  (XO(I)  ,  1=1, Nl) 

INT29550 

c 

WRITE  (6  ,  298) 

INT29560 

c 

WRITE  (6  ,  300)  ( Y0( I)  ,  1=1, Nl) 

INT29570 

c 

INT29580 

RETURN 

INT29590 

c  -  ■ 

INT29600 

100 

FORMAT(7I5) 

INT29610 

200 

F0RMAT(6F10.  0) 

INT29620 

290 

FORMAT ( / ,  ’  ORIGINAL  COORDINATES’ '  X/C’) 

INT29630 

295 

FORMATC/,'  INTERPLATED  COORDINATES’ ,/, '  X/C’) 

INT29640 

298 

FORMAT ( ’  Y/C') 

INT29650 

300 

FORMAT( 6F10.  6) 

INT29660 

C  -  • 

INT29670 

END 

INT29680 

c 

INT29690 

SUBROUTINE  STAGR( N , STAG , XO , YO , XSTGR , YSTGR) 

INT29700 

c 

INT29710 

DIMENSION  X0( 100) ,Y0( 100) ,XSTGR( 100) ,YSTGR( 100) ,DS( 100) 

INT29720 

c 

INT29730 

XOTE  =  0.5  *  (XO( l)+XO(N) ) 

INT29740 

YOTE  =  0.  5  *  ( YO( 1 )+YO(N) ) 

INT29750 

DS( 1)  =  SQRT( (X0( 1) -X0TE)**2  +  ( Y0( 1) -Y0TE)**2) 

INT29760 

DSM  =  DSC  1) 

INT29770 

DO  10  I  =  2  ,  N 

INT29780 

DS( I )  =  SQRT( C  X0( I ) -XOTE) **2  +  (Y0(I)-Y0TE)**2) 

INT29790 

IF  C DS( I )  ■  LT.  DSM)  GOTO  10 

INT29800 

IM  =  I 

INT29810 

DSM  =  DS( I) 

INT29820 

10 

CONTINUE 

INT29830 

c 

INT29840 

YYY  =  YOTE-YOCIM) 

INT29850 

XXX  =  XOTE-XO(IM) 

INT29860 

IF  C  YYY  .EQ.  0.0  .AND.  XXX  .  EQ.  0.0)  THEN 

INT29870 

ANG  =  0.  0 

INT29880 

ELSE 

INT29890 

ANG  =  ATAN2C Y"  I  ,XXX) 

INT29900 

END  IF 

INT29910 

ANG  =  ANG  +  STAG 

INT29920 

c 

INT29930 

COS AN  =  COS (ANG) 

INT29940 

SI NAN  =  S INC ANG) 

INT29950 

DO  20  I  =  1  ,  N 

INT29960 

c 

YY  =  YOCI)-YO(IM) 

INT29970 

c 

XX  =  XOCI)-XO(IM) 

INT29980 

c 

IF  (YY  .EQ.  0.0  .AND.  XX  .  EQ.  0.0)  THEN 

INT29990 

c 

ANGCO  =  0.0 

INT30000 

c 

ELSE 

INT30010 

c 

ANGCO  =  ATAN2(YY,XX) 

INT30020 

c 

END  IF 

INT30030 

XSTGR( I )=  XO( I )*COSAN  +  YO(I)*SINAN 

INT30040 

YSTGR( I )=  YO(I)*COSAN  -  XO(I)*SINAN 

INT30050 

20 

CONTINUE 

INT30060 
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RETURN' 

END 
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APPENDIX  B.  C4  CASCADE 


A.  EXPERIMENTAL  RESULTS 

The  experimental  results  of  the  C4  cascade  were  obtained  directly  from  professor 
GJ.  Walker,  University  of  Tasmania,  Tasmania,  Australia,  who  performed  these  exper¬ 
iments. 

The  results  of  the  boundary  layer  measurements  of  the  C4  cascade  are  given  below 
at  four  inlet  angles:  34.1°,  36.3°,  45.6°,  and  47.7°.  The  Reynold  numbers,  based  on  the 
chord  and  the  upstream  velocity,  are  200000,  191000,  173000  and  171000  respectively. 
The  results  given  in  the  following  tables  include  the  displacement  thickness  (<5*),  the 
shape  factor  { H )  and  the  local  free  stream  velocity  (L'E). 


Table  1.  EXPERIMENTAL  RESULTS  AT  INLET  ANGLE  OF  34.1° 


X  c 

S'  [10~3  FT] 

H 

L'E  [FT/SEC] 

0.4 

4.9 

2.48 

16S.37 

0.5 

6.28 

2.61 

167.35 

0.6 

S.79 

3.24 

158.31 

0.7 

10.83 

3.63 

149.27 

o.s 

16.63 

3.79 

147.13 

0.9 

16.19 

1.89 

143.79 

Table  2.  EXPERIMENTAL  RESULTS  AT  INLET  ANGLE  OF  36.3° 


x  c 

cf[I0-a  FT] 

H 

UE  [FT  SEC] 

0.4 

5.43 

2.55 

161.63 

0.5 

7.09 

2.70 

157.59 

0.6 

10.3 

3.34 

148.78 

0.7 

12.63 

3.78 

139.87 

0.8 

14.84 

2.78 

135.01 

0.9 

16.43 

1.76 

133.23 
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Table  3.  EXPERIMENTAL  RESULTS  AT  INLET  ANGLE  OF  45.6° 


\  c 

5*  [  1 0-3  FT 

II 

UE  [FT  SEC] 

0.4 

s.os 

2.5S 

137.SS 

0.5 

9. S3 

2.41 

133.70 

0.6 

12.35 

2.33 

122. IS 

0.7 

12.9S 

1.97 

114.93 

O.S 

19.44 

1.90 

111.77 

0.9 

27.69 

1.92 

109.26 

Table  4.  EXPERIMENTAL  RESULTS  AT  INLET  ANGLE  OF  47.7° 


x  c 

cMlO-4  FT] 

H 

UE  [FT  SEC] 

0.4 

S.87 

2.24 

130.64 

0.5 

10.27 

2.19 

124.76 

0.6 

14.31 

2.08 

116.86 

0.7 

16.45 

1.S7 

106.72 

0.8 

24.16 

1.S2 

103.75 

0.9 

36.00 

2.01 

102. IS 

l 
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The  results  of  the  measurements  of  the  velocity  profiles  in  the  boundary  layer  at  two 
inlet  angles.  34.1°  and  36.3'  at  5<)"0  chord  are  given  below. 


Table  5.  VELOCITY  PROFILES  AT  50%  CHORD. 


y 

P  =  36.3° 

p  =  34.1° 

0.0 

0.0 

0.0 

2.3 

0.172 

0.20S 

3.7 

0.270 

0.327 

6.2 

0.469 

0.534 

8.6 

0.666 

0.728 

1 1 .0 

0.794 

0.867 

13.4 

0.S91 

0.933 

IS. 3 

0.9S2 

0.985 

23.2 

1 .000 

1.000 

B.  C4  CASCADE  COORDINATES 


DIMENSION  X(0:  100) ,XU(0:  100) ,XL(0:  100) ,YU(0:  100) ,YL(0: 100) 
DATA  A1 ,A2 , A3 ,A4/0.  15492,0. 06563,0. 2528,0. 2811/ 

DATA  31 , B2 , B3 ,B4/0. 03866,0.  07871,0. 1467,0. 03448/ 

PI  =  ACQSC-l. 0) 

C  READ  (5,800)  NMAX 

800  FORMAT  (15) 

NMAX=33 

C  READ  (5,810)  (X( I) ,1=0, NMAX) 

810  FORMAT  (6F10. 6) 

DO  50  1=0, NMAX 

X(I)  =  (1.  0-COS(PI*I/NMAX))/2. 

50  CONTINUE 

DO  100  1=0, NMAX 

SRT  =  SQRT((0.  5/SIN(PI/12))**2-(0. 5-X(I))**2) 

YC  =  -0. 5/TAN(PI/12)  +  SRT 
DY  =  ATAN((0.  5-X( I) )/SRT) 

IF  (X(I).  LT.  0.  3)  THEN 

YT  =  A1*SQRT(X( I ) )  -  A2*X(I)  -  A3*X(I)**2  +  A4*X(I)**3 
ELSE 

YT  =  B1  +  B2*X( I)  -  B3*X(I)**2  +  B4*X(I)**3 
END  IF 


C4  00010 
C4  00020 
C4  00030 
C4  00040 
C4  00050 
C4  00060 
C4  00070 
C4  00080 
C4  00090 
C4  00100 
C4  00110 
C4  00120 
C4  00130 
C4  00140 
C4  00150 
C4  00160 
C4  00170 
C4  00180 
C4  00190 
C4  00200 
C4  00210 
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o  o 


YU(I)  =  YC  +  COS(DY)*YT  C4  00220 

YL( I )  =  YC  -  COS(DY)*YT  C4  00230 

XU(I)  =  X(I)  -  SIN(DY)*YT  C4  00240 

XL( I )  =  X(I)  +  SIN( DY)*YT  C4  00250 

100  CONTINUE  C4  00260 

WRITE  (6,900)  ( I ,X( I ) ,XU( I ) , YU( I ) ,XL( I ) , YL( I ) , 1=0 , NMAX)  C4  00270 

900  FORMAT  ( 15 ,4X, F10. 6 ,4X , 2FI0. 6 ,4X , 2F10. 6)  C4  00280 

WRITE  (1,910)  ( XL( I ) , I=NMAX , 0 , - 1 ) , ( XU( I ) , 1=1 , NMAX)  C4  00290 

WRITE  (1,910)  (YL(I),I=NMAX,0,-1),(YU(I), 1=1, NMAX)  C4  00300 

910  FORMAT  (6F10. 6)  C4  00310 

STOP  C4  00320 

END  C4  00330 


ro 
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